Tools for quantum network design

Quantum networks will enable the implementation of communication tasks with qualitative advantages with respect to the communication networks we know today. While it is expected that the first demonstrations of small scale quantum networks will take place in the near term, many challenges remain to scale them. To compare different solutions, optimize over parameter space and inform experiments, it is necessary to evaluate the performance of concrete quantum network scenarios. Here, we review the state of the art of tools for evaluating the performance of quantum networks. We present them from three different angles: information-theoretic benchmarks, analytical tools, and simulation.

IV. Simulation tools 24 A. Simulation methodology 24 B. The challenge of simulating quantum networks 25 C. Existing simulation platforms 25

V. Outlook 26
Data availability statement 26

I. INTRODUCTION
Communication networks are the basis of our connected society. In particular, the internet that we know today is a conglomerate of physical links of very different characteristics. Some are highly reliable and persistent, while some are very noisy and dynamic. From this basic resource, the network builds services, such as reliable remote transmission, multicast, and many others.
Network engineering involves challenges of very different nature. At the hardware level, there might be the request of replacing a device element or of scaling an architecture from tens to thousands of nodes. At the network level, there might be the need to develop a new protocol. At the application level, new applications with different traffic patterns and requirements on the network might be deployed. To address these challenges and be able to take appropriate measures, it is necessary to evaluate solutions and then compare against alternatives and benchmarks in a quantitative way.
Performance can be evaluated roughly in three different ways: experimentally, by means of simulation, and analytically. Experimental methods evaluate the performance with hardware. They can range from prototyping to fully-fledged field test deployments. The advantage of experimental methods is that they give the most accurate answer to a concrete scenario, but they have several drawbacks. Notably, they can be extremely costly and offer little flexibility to change parameters. Moreover, they might be unavailable to evaluate directions of future research, i.e. if the hardware does not yet exist. Simulation and analytical methods require models of the different network elements. A first observation is that the quality of the evaluation depends on the accuracy with which the models capture the behavior of each element. The difference between simulation and analytical methods is that simulation methods imitate the behavior of the individual elements in the network and their interactions, while analytical methods compose analytical models representing some characteristic of interest to determine the aggregate performance without actually replicating the actions of the elements. Both simulation and analysis can tackle a broad range of scenarios at the cost of a less accurate evaluation.
The choice of the evaluation method and the concrete tool depend on many factors. Some of them are cost, flexibility, accuracy, and validation.
Moreover, the methods are not exclusive. For instance, a streamlined model for mathematical analysis can be validated by simulation or experiment. Or alternatively, one can use an analytical model to validate a new simulation approach. In general, cross-validation is a well known method for increasing the credibility and reliability of an evaluation tool 1 .
In principle, quantum networks enable the implementation of tasks beyond the reach of classical networks such as key distribution 2,3 , clock synchronization 4 , increasing the baseline of telescopes 5 , and many others 6 . However, while they share similar high-level challenges with their classical counterpart, there are noticeable differences. On the one hand, quantum information can not be copied, reducing the applicability of classical approaches for long-distance communications. On the other hand, entanglement between two parties is only useful if it is clear which particles at one location are entangled with which ones at the distant location. In turn, this places novel constraints to the network architecture.
Motivated by recent experimental advances promising the deployment of the first quantum networks, the past years have seen a wealth of research dedicated to analytical and simulation tools for evaluating the performance of quantum networks. Our goal in this paper is to review them and make a unified description of all existing resources. We undertake this task in three steps.
First, in Section II, we review analytical information-theoretic tools. Given a model for communication channels and some precise definition of the available resources, information theory allows us to bound the optimal rate at which a communications task can be performed. This abstract optimal rate, called the channel capacity, typically can only be achieved in a very ideal setting, for instance as the number of channel uses tends to infinity. Hence, information-theoretical bounds represent a fundamental performance limit and can be used to benchmark concrete implementations. One important example of their role is the evaluation of quantum repeater implementations. Quantum repeaters enable, in principle, the transmission of quantum information over arbitrarily long distances with non-vanishing rates. While there has been strong experimental progress, it is unclear what performance is necessary to label an implementation as a quantum repeater. Information theory provides a clear tool for evaluating implementations. If two parties are connected by a quantum channel, their communication rate is bounded by the channel capacity. If the two parties place a device between them and, with the help of this device, they achieve a communication rate above capacity this certifies that the device is qualitatively improving the communication capabilities of the two parties. Second, in Section III, we review analytical tools for evaluating the performance of concrete protocols. We focus our analysis in protocols for distributing entanglement as, once the entanglement has been distributed, it can be consumed for implementing different tasks. Given a protocol, the key performance parameters are the time it takes to distribute the entanglement and the quality of the entanglement. The main difficulty for this analysis is that many protocols rely on probabilistic elements and these elements can be stacked in nontrivial ways to construct the entanglement distribution protocol. In consequence, both the waiting time and the entanglement quality are random variables. In spite of this, as long as the physical models are simple enough, it is possible to approximate, and sometimes calculate, relevant parameters of these random variables for a large class of protocols including entanglement swapping, distillation, and cut-offs.
Meanwhile, simulation is the tool of choice for understanding different types of complex behavior in classical networks. Until recently, there existed very few options for quantum network simulations. In Section IV, we discuss the role of simulation and the existing platforms. For this, first, we outline the methodology for evaluating the network performance based on simulation and then discuss the particularities of quantum networks together with a selection of existing tools for quantum networks.
We end the review with a summary of the tools presented and an outlook in Section V.
Let us clarify the scope of the review. We review tools for evaluating the performance of quantum networks as the key element that will enable the design of quantum networks. We leave the actual design of quantum networks [7][8][9][10][11][12][13][14][15][16][17] beyond the scope of this review. We also omit the closely related topics of quantum networks for creating correlations [18][19][20][21][22] and of complex quantum networks 23,24 . Finally, there is a large literature regarding communications over direct links connecting two adjacent nodes. The fundamental limits over such a "network" for a task are given by the associated channel capacities, while the performance of a concrete protocol can be estimated in many different ways. There exist a number of references reviewing the fundamental and practical methods to evaluate a direct connection [25][26][27] . For this reason, we have chosen to restrict the scope to tools that are particular to networks beyond direct connections.

II. FUNDAMENTAL LIMITS FOR ENTANGLEMENT DISTRIBUTION OVER NETWORKS
Even in the classical case, physical channels connecting two end points can have a very complicated behavior. However, for practical purposes, many communication channels can be abstracted into a stochastic map that takes a symbol at the input and places a symbol at the output following some probability distribution. One can then study the usefulness of this channel for a communication task. The usefulness is measured by the number of times that the task can be achieved divided by the number of uses of the channel. The optimal ratio maximized over all possible communication protocols without computational restrictions is called the capacity of the channel. While this is a very abstract notion, it is useful for benchmarking practical communication schemes. Moreover, the capacity of a classical channel can be easily approximated numerically 28,29 . Quantum channels can also be abstracted by a mathematical object (see below II A), and the capacity of a quantum channel for different communication tasks can be defined. However, it is presently unknown how to compute the capacity for a general channel and in many cases only bounds are available. Fortunately, for some very relevant channels such as the lossy bosonic channel, the capacity is known 30 . Building on the tools for point-to-point channels, a number of recent results allow to also bound the capacities for many communication tasks over networks, even including multi-user scenarios. In the rest of the section, we present a general description of such tasks, upper/lower bounds on their capacities, and efficient ways to evaluate the bounds.

A. Networks as graphs
A quantum network is composed of quantum information processing nodes and quantum channels. A quantum information processing node may represent a client, a quantum computer, or a quantum repeater node, but, in theory, it is regarded as an abstract node which can perform arbitrary local operations (LO) allowed by quantum mechanics 27,31 . This local operation is, in general, a stochastic process, described by a completely-positive (CP) map. More precisely, a local operation at node X receives a quantum statê ρ as an input, and returns a quantum stateσ k with probability p k as an output, whereσ k :=M X kρ (M X k ) † /p k with Kraus operators {M X k } k satisfying ∑ k (M X k ) †MX k =1 X and p k = Tr[(M X k ) †MX kρ ]. On the other hand, a quantum channel is a device to convey a quantum system from one node to another node, such as an optical fibre, a superconducting microwave transmission line, or an optical free space link. In theory, any quantum channel N X→Y from a node X to a node Y is described by a completely-positive trace-preserving (CPTP) map. In particular, quantum channel N X→Y for an input stateρ is described by N X→Y (ρ) := Tr E [Û XE→Y E (ρ ⊗ |0 0| E )(Û XE→Y E ) † ] with a unitary operatorÛ XE→Y E on Hilbert space H X ⊗ H E = H Y ⊗ H E and a state |0 E of auxiliary system E.
The structure of a quantum network, which might be a subnetwork of a larger quantum network in general, can be specified in terms of a graph (Fig. 1a). A quantum information processing node in a quantum network is regarded as a vertex v in a graph, while a quantum channel N v→v to send a quantum system from a node v to another node v is associated with a directed edge e = v → v (or, simply, e = vv , where the former v and latter v are the tail t(e) and head h(e) of the edge e, respectively) in the graph. With this rule, a given quantum network is specified by a graph G = (V, E) with a set V of vertices and a set E of directed edges, in such a way that vertices in the set V and directed edges in the set E correspond to quantum information processing nodes and quantum channels in the quantum network, respectively.
In addition to a quantum network, we may use a classical communication network, such as a telephone network or the Internet. In general, this classical communication network may have a restriction, for instance, such that a quantum information processing node is not connected to another node with a classical communication channel, or such that classical communication between quantum information processing nodes is restricted to only one way. In practice, such a restriction about classical communication should be taken care of. However, if we are interested in the ultimate performance of a quantum network, or if classical communication is supposed to be much cheaper than quantum communication in the era of quantum networks, then it may be reasonable to assume that classical communication between any nodes is free. Since LO are also free within every quantum information processing node, this leads to an assumption that local operations and classical communication (LOCC) among all quantum information processing nodes are free.

B. Communication tasks
A major goal of the quantum internet is to distribute entanglement among clients in a quantum network. This is because once clients share an entangled state, a communication task can be executed by consuming the entanglement, under LOCC, that is, without using the quantum network anymore. What kind of entanglement is required depends on the task the clients want to perform. For example, if a client A(∈ V ) wants to send an unknown quantum state with dimension d to another client B(∈ V ) via quantum teleportation 33 , they need to share a bipartite maximally entangled state |Φ d AB : Even if the client A wants to send bits to the client B secretly, since the maximally entangled state provides them log 2 d bits of secret key for the one-time pad, sharing the maximally entangled . . ,C 9 } is regarded as a quantum network associated with a graph G = (V, E), where edge e ∈ E is used to specify the existence of a quantum channel N e in the network. A cut ∂ (V A ) is composed of green edges that separate the network into sets V A = {A,C 1 ,C 2 , . . . ,C 5 }, shaded in red, and . . ,C 9 }, shaded in green. (b) Schematic description of the protocol performed at a physical layer of a linear network. In the first round channel N e 1 is used, followed by LOCC, providing output k 1 . In the second round, depending on outcome k 1 in the previous round, channel N e k 1 is used followed by LOCC providing output k 2 and so on. Reprinted with permission from K. Azuma . More generally, one could think of a group of clients C := C 1 C 2 · · ·C n (⊂ V ) (C 1 ,C 2 , . . . ,C n ∈ V ) who want to share a multipartite entangled state, such as a Greenberger-Horne-Zeilinger (GHZ) state 36 which can be used for secret sharing 37 , a clock synchronization protocol 4 as well as quantum conference key agreement 38 . In the case of quantum conference key agreement there is, again, a larger class of so-called multipartite privatesγ C d , from which the group of clients can obtain log d bits of conference key 38 . Such states are of the formγ C d : is an arbitrary controlled unitary operation in the formÛ C 1 C 2 ···C n twist := Depending on the communication task that clients want to perform, clients C may ask a provider to supply an resource entangled stateτ C d for the task. For example,τ C d may be a bipartite maximally entangled state |Φ d C , a private stateγ AB d , a GHZ state |GHZ d C or a multipartite private stateγ C d . The amount of the entangled resource present inτ C d is quantified by the parameter d. Namely, log 2 d is the number of entangled qubits in a maximally entangled or GHZ state, or the number of secret bits available in a bi-or multipartite private state. The case of cluster states and general graph states will be discussed in Section II E 2. According to the request from the clients C, the provider prepares an available quantum network, associated with a graph G = (V, E) such that C ⊂ V . This is because entanglement cannot be generated only with LOCC in general (or, precisely, entanglement is defined 39 as a state which cannot be generated by LOCC, that is, which cannot be described by the form of fully separable states ∑ i q i v∈Vρ v i , where {q i } i is a probability distribution andρ v i is a quantum state of node v ∈ V ), and thus, the provider needs to use a quantum network. Without loss of generality, a protocol of the provider which provides a final stateρ V k l close to a target stateτ C d k l with probability p k l is described as follows (see also Fig. 1b). A provider determines an available quantum network, associated with a graph G = (V, E) such that C ⊂ V . The protocol starts by preparing the whole system in a fully separable stateρ V 1 and then by using a quantum channel N e 1 with e 1 ∈ E. This is followed by arbitrary LOCC among all the nodes, which gives an outcome k 1 and a quantum stateρ V k 1 with probability p k 1 (because LOCC is based on local CP maps and sharing their outcomes by classical communication). In the second round, depending on the outcome k 1 , a node may use a quantum channel N e k 1 with e k 1 ∈ E, followed by LOCC among all the nodes. This LOCC gives an outcome k 2 and a quantum stateρ V k 2 k 1 with probability p k 2 |k 1 . Similarly, in the i-th round, according to the previous outcomes k i−1 := k i−1 . . . k 2 k 1 (with k 0 := 1), the protocol may use a quantum channel N e k i−1 with e k i−1 ∈ E, followed by LOCC providing a quantum statê ρ V k i with a new outcome k i with probability p k i |k i−1 . Eventually, after a finite number of rounds, say after an l-th round, the protocol presents clients C withρ 1 ≤ ε for ε > 0, with probability p k l := p k l |k l−1 · · · p k 3 |k 2 p k 2 |k 1 p k 1 , where X 1 is the trace norm defined by X 1 := Tr √X †X . Important quantities related to the evaluation of a protocol are (1) how many times quantum channels N e are used in the protocol, (2) how valuable the stateτ C d given to clients C is, and (3) how much error ε is allowed. As for (1), given that a general protocol could work in a probabilistic manner, l e := ∑ l−1 i=0 δ e,e k i k i = ∑ l−1 i=0 ∑ k i p k i δ e,e k i with the Kronecker delta δ i, j andl E := ∑ e∈El e are important quantities, becausel e represents the average number of times quantum channel N e is used in the protocol andl E represents the average number of total channel uses. Thus, {l e } e∈E can be a set to characterize the protocol. If we introduce an average usage frequencȳ f e :=l e /l E (≥ 0) satisfying ∑ e∈Ef e = 1, a set {l E , {f e } e∈E } can be used to characterize the protocol. Or, by introducing an average usage rater e :=l e /l(≥ 0), a set {l, {r e } e∈E } can also be used to characterize the protocol, where l is a parameter like time (because it has no explicit relation withl e in contrast tol E ). As for (2), in general, it is not so obvious to define a quantity for evaluating the value of the target stateτ C d k l , because the value depends on the final application in general and may be determined by a complex function of the target stateτ C d or its parameter d. However, if the target stateτ C d is GHZ state |GHZ d C or private state γ C d , the value of the target state can be quantified by log 2 d showing how many qubits or private bits can be sent by consuming the target state under LOCC. As for (3), error ε in terms of the trace distance quantifies how far from a target stateτ C k l the final stateρ C k l is. In general, this trace distance ρ −σ 1 is related 40 with another well known measure, the fidelity F(ρ,σ ) := ρ √σ 1 , as Although the fidelity is a useful and intuitive tool, the trace distance is used often in theory for quantum communication, because the trace distance has universal composability 41 , in contrast to the fidelity.

For two-party communication
The upper bounds on how efficiently maximally entangled states and private states can be distributed for two clients have been derived by Pirandola 42,43  follows.
An entanglement measure 39 E between parties X and Y does not increase on average under any LOCC between them and is zero for any separable state. In addition to this normal property as an entanglement measure, suppose that the entanglement measure E satisfies the following two properties: is called the entanglement of channel 30,44,46,47 . With an entanglement measure E satisfying Properties (i) and (ii), an upper bound is given in the following form: For any protocol which provides clients A(∈ V ) and B(∈ V ) with maximally entangled state |Φ d k l AB (or private stateγ AB d k l ) with probability p k l and error ε > 0, by using a quantum network associated with graph G = (V, E), it holds is the set of edges which connect a node in V A and a node in V \ V A (Fig. 1). Since Eq. (1) holds for any choice of V A , the minimization of the right-hand side of Eq. (1) over the choice of V A gives the tightest bound on the left-hand side. In general, there might be cases where we have no computable expressions for pointto-point capacities, the best characterization which is simply the optimal rate achievable with a protocol over any number of channel uses. In consequence, it is unclear how to proceed to optimize additionally over the average usel e of each channel N e . In contrast, the upper bound (1) is always estimable, once we are given {E (N e )} e∈E or their upper bounds, because it depends only on a single use of channel N e , rather than multiple uses of it.
Examples of entanglement measure E satisfying Property (i) 45 for any stateρ A , where d = dim H A . Examples of these channels are ones which "commute" with the unitary correction in the quantum teleportation. In particular, the quantum teleportation 33 from system A in an unknown stateρ A to system A is implemented by generalized Bell measurement for any statê σ A and any outcome i, then the channel N A→B is teleportation simulable. For example, depolarizing channels, phaseflip channels, and lossy bosonic channels are teleportation simulable, while amplitude damping channels are not teleportation simulable.
In contrast to cases of point-to-point quantum communication (that is, cases corresponding to a graph G with V = {A, B} and E = {A → B}), the definitions of quantum/private capacities in a quantum network themselves may have variety 63,64 .
A rather general definition is for given {q e } e∈E with q e ≥ 0, where n ≥ 0, Λ(n, {q e } e∈E , ε) is a set of protocols which provide clients A(∈ V ) and B(∈ V ) with maximally entangled state |Φ d k l AB (or private statê ) with probability p k l and error ε > 0 by using quantum channels N e nq e times on average at most (i.e.,l e ≤ nq e ), and C represents a network quantum capacity Q (private capacity P) per time. Notice that any protocol with lr e ≤ nq e belongs to the set Λ(n, {q e } e∈E , ε) of protocols. We can capture alternative definitions of capacity by putting additional constraints on {q e } e∈E . Let us introduce two particularly relevant alternative capacity definitions.
First, we can consider a scenario where the users can optimize over the usage frequency of individual channels and define the capacity as the optimal rate normalized by the total number of channel uses: where C represents a network quantum capacity Q (private capacity P) per the average total number of channel uses and we have added the additional constraint: ∑ e∈E q e = 1. Notice that any protocol withl Ef e ≤ nq e belongs to the set Λ(n, {q e } e∈E , ε) of protocols with ∑ e∈E q e = 1. Second, we can consider the usage of the whole network as the unit resource. This network use metric can be captured by setting q e = 1 for all e ∈ E. This quantity is equivalent to the one introduced by Pirandola 43 for multi-path routing protocols with a so-called flooding strategy.
Finally, we note that Pirandola also introduced the single path per network use capacity 43 . This quantity is not known to correspond to a restriction on {q e } e∈E . For additional details on this quantity we refer to the original reference 43 and to a reference 64 for a comparison with the channel use metric.
For all of these capacities except for the single path per network use, Eq. (1) directly gives upper bounds, for C = Q, P. Also note that holds, because the maximally entangled state |Φ d AB is a special example of the private statesγ d AB . An upper bound in the form (1) was originally derived for point-to-point quantum communication, that is, for cases corresponding to a graph G with V = {A, B} and E = {A → B}. The first useful upper bound on the private/quantum capacity is derived by Takeoka et al. 44 , with the use of the squashed entanglement E sq , which applies to any point-to-point quantum channel N A→B . This TGW bound has no scaling gap with secure key rates (per pulse) of existing point-to-point optical QKD protocols, showing that there remains not much room to improve further the protocols in terms of performance. An alternative bound is derived by Pirandola et al. 30 , with the use of the relative entropy E R of entanglement, to obtain tighter bounds for teleportation simulable channels. Indeed, this PLOB bound coincides with the performance of achievable point-to-point entanglement generation protocols over erasure channels, dephasing channels, bosonic quantum amplifier channels, and lossy bosonic channels, implying that it represents the quantum/private capacities, i.e., for those channels. For instance, in the case where the channel N A→B is a single-mode lossy bosonic channel with a transmittance η (0 ≤ η ≤ 1), However, the PLOB bound can be applied only to teleportation simulable channels, in contrast to the TGW bound. As a general upper bound applicable to arbitrary channels like the TGW bound, Christandl and Müller-Hermes derive an upper bound 46 , with the use of the max relative entropy E max of entanglement, which is also shown to be a strong converse bound, i.e. the error tends to one exponentially once the rate exceeds the bound. Further detail on the comparison between the upper bounds can be found in the review paper by Pirandola et al. 65 . Then, those upper bounds for point-to-point quantum communication were lifted up to ones for general quantum networks associated with a graph G = (V, E), by Pirandola 42,43 , Azuma et al. 32 and Rigovacca et al. 45 . The main idea here is to regard a quantum internet protocol as a point-to-point quantum communication protocol between a party having the full control over V A (⊂ V ) and a party having the full control over V \V A . In particular, Pirandola derives 42,43 an upper bound (1) with E = E R for any quantum network composed of teleportation simulable channels, with the use of the relative entropy E R of entanglement, by following the PLOB bound 30 , while Azuma et al. provide 32 an upper bound (1) with E = E sq for any quantum network composed of arbitrary quantum channels, with the use of the squashed entanglement E sq , by following the TGW bound 44 . To make Pirandola's relative entropy bound applicable to arbitrary quantum networks like the upper bound of Azuma et al., Rigovacca et al. present 45 a 'hybrid' relative entropy bound for any quantum network composed of arbitrary quantum channels, by using an inequality introduced by Christandl and Müller-Hermes 46 . This bound corresponds to the case where E (N e ) in Eq. (1) is assumed to be where N e ∈ S indicates that quantum channel N e is teleportation simulable. Pirandola gives 42,43 an achievable protocol that works over a quantum network composed of teleportation simulable channels, while Azuma et al. present 63 an abstract but general achievable protocol working over any quantum network. Here we focus on the latter protocol. In particular, Azuma et al. give an achievable protocol over a general quantum network associated with a graph G = (V, E), by assuming that we have an achievable entanglement generation protocol over every channel N e . In particular, they assume that there is an entanglement generation protocol (including entanglement purification) over quantum channel N e which provides a stateρ e δ -close to l e R e δ copies of a qubit maximally entangled state |Φ 2 AB , called a Bell pair, by using quantum channel (N e ) ⊗l e , that is, ρ e − (|Φ 2 Φ 2 | e ) ⊗l e R e δ 1 ≤ δ with δ ≥ 0 (for large l e ). If this entanglement generation protocol runs over every channel N e , the network can share the state If we regard each Bell state |Φ 2 Φ 2 | e of the Bell-pair network e∈E (|Φ 2 Φ 2 | e ) ⊗l e R e δ as an undirected edge e with the same two ends of e, we can make a multigraph G = (V, E ) with the set of E of such undirected edges e . In this multigraph G , if we find a single path between nodes A and B, called AB-path, then we can give a Bell pair to the nodes A and B by performing the entanglement swapping at the vertices along the AB-path. Hence, the number of Bell pairs which are provided to the nodes AB equals to the number of edge-disjoint AB-paths. Menger's theorem in graph theory 66,67 tells us that the number M of edge-disjoint AB paths in a multigraph is equivalent to the minimum number of edges in an AB cut (see Section II C 1 for this theorem). Thus, M = min V A ∑ e∈∂ (V A ) l e R e δ for the multigraph G . If we perform the entanglement swapping along M edge-disjoint AB paths, we have from Eq. (11). Therefore, this protocol, called the aggregated quantum repeater protocol, supplies nodes AB with min V A ∑ e∈∂ (V A ) l e R e δ Bell pairs with error |E|δ . For fixed q e = l e /l(> 0) (or q e = l e /l E (> 0) with l E = ∑ e∈E l e ), if the assumed entanglement generation protocol over every channel N e can achieve quantum capacity, that is, δ → 0 and R δ → Q(N e ), by regarding l (l E ) as n and then taking n → ∞, the aggregated quantum repeater protocol gives lower bounds on the network quantum capacities of protocols Λ(n, {q e } e∈E , ε) as max If a quantum network is composed only of quantum channels N e that satisfy Eq. (9), as Pirandola has considered 42,43 , lower bounds in Eqs. (13) and (14) equal to upper bounds in Eqs. (5) and (6) with E = E R , respectively. For instance, in the case of a quantum network composed only of lossy bosonic channels, the aggregated quantum repeater protocol achieves the private/quantum capacities of the network 63 .

For multiple-pair communication
In general, a provider may be asked to concurrently distribute multiple pairs of bipartite maximally entangled states or bipartite private states, i.e.,τ C d : (1) , d (2) , · · · ). The upper bounds for this problem are given 42,68 by Pirandola for any quantum network composed of teleportation simulable channels, with the use of the relative entropy E R of entanglement, while they are given by Bäuml et al. 69 is the set of edges which connect a node in V and a node in V \V , and I(V A ) is the set of indices i whose corresponding pair Capacities for multiple-pair communication could have more variety than those for two-party communication, because there are many pairs of clients 64

. A possible definition is
) with probability p k l and error ε > 0 (i.e., , by using quantum channels N e nq e times on average at most (i.e.,l e ≤ nq e ), and C worst represents a worst-case network quantum capacity Q worst (private capacity P worst ) per time. Notice that any protocol with lr e ≤ nq e belongs to the set Λ(n, {q e } e∈E , ε) of protocols. The use of objective function min i { log 2 d (i) k l k l /n} i means that a protocol is evaluated with the least achievable rate that is guaranteed for any client pair. By putting an additional constraint ∑ e∈E q e = 1 on the set {q e } e∈E , the capacity per time is transformed into one per the average total number of channel uses: . (17) Notice that any protocol withl Ef e ≤ nq e belongs to the set Λ(n, {q e } e∈E , ε) of protocols with ∑ e∈E q e = 1. Another possible definition is for given {q e } e∈E with q e ≥ 0 and given {s i } i with s i ≥ 0, where C weight represents a weighted multi-pair network quantum capacity Q weight (private capacity P weight ) per time. The use of objective function ∑ i s i log 2 d (i) k l k l /n means that a protocol is evaluated with the average of rates for client pairs A i B i . By putting an additional constraint ∑ e∈E q e = 1 on the set {q e } e∈E , another possible definition is where C weight represents a weighted multi-pair network quantum capacity Q weight (private capacity P weight ) per the average total number of channel uses. As an important special case of C weight , we may also consider total throughput C total defined as (20) and Eq. (15) gives upper bounds on the capacities. Let R i := log 2 d (i) k l k l /n. If we multiply Eq. (15) by n −1 with n ≥ 0 and take the limit of n → ∞ and ε → 0, it gives upper bounds on R i for every i, on R i + R j for every i = j, on R i + R j + R k for every different i, j, k, · · · , and on ∑ i R i , by changing the choice of V . For given {q e } e∈E , if we minimize every upper bound over possible choices of V for it, the minimized upper bounds give a minimal polytope in the Cartesian coordinates R := (R 1 , R 2 , · · · ) which includes all possibly achievable points of vector R. Since objective functions in the definitions of , we need to start again from finding the minimal polytope and then to maximize the objective function. Instead of these recipes including not only minimization but also maximization, Bäuml et al. provide 64 a computationally efficient method to obtain upper bounds, which are explained in Section II D.
To find lower bounds on the capacities, we need to construct a protocol. A simple way to make such a protocol is to invoke the idea of the aggregated quantum repeater protocol 63 . In particular, by running entanglement generation over every quantum channel N e , the network shares a quantum state e∈Eρ e ε-close to e∈E (|Φ 2 Φ 2 | e ) ⊗l e R e δ associated with an undirected multigraph G = (V, E ), where each undirected edge e ∈ E with the same two ends of e corresponds to a Bell pair |Φ 2 e between the ends of e . Then, if we find a path between A i and B i in the undirected multigraph G , we can provide a state close to a Bell pair |Φ 2 A i B i by performing entanglement swapping along the A i B i -path. Thus, depending on an objective function in a definition of a capacity, by finding out edge-disjoint paths between pairs A i B i in the undirected mutligraph G properly, we can give a state close Along this recipe, Bäuml et al. provide 64 a computationally efficient way to obtain lower bounds on the capacities, which are explained in Section II D.

For multipartite communication
The upper bound 69  Indeed, they have considered a quantum broadcast network which might include quantum broadcast channels as well. A quantum broadcast channel is a device to distribute a quantum system from a single node to multiple nodes, and thus, the edge e of a quantum broadcast channel N e is a directed hyperedge which has a single tail t(e) ∈ V and multiple heads h(e)(⊂ V ). This implies that a quantum broadcast network is associated with a directed hypergraph G = (V, E) with a set V of vertices and a set E of directed hyperedges. Besides, the goal of a quantum internet protocol working over such a quantum broadcast network is to distribute multipartite GHZ statesτ C d : Then, the upper bound is described as follows, with the use of multipartite squashed entanglement E P sq : Any protocol which provides every set C i of clients with GHZ state |GHZ probability p k l and error ε > 0, by using a quantum broadcast network associated with a directed hypergraph G = (V, E), follows for any partition P = P 1 : P 2 : · · · : P k of the set V , where is the set of hyperedges whose tail t(e) and heads h(e) all belong to one set of P, n C i |P = 0 for the case of |{l ∈ {1, . . . , k} : P l ∩C i = / 0}| < 2 and n C i |P = |{l ∈ {1, 2, . . . , k} : P l ∩ C i = / 0}| for the case of |{l ∈ {1, 2, . . . , k} : The multipartite squashed entanglement of broadcast channel N e is defined by E P sq (N e ) := max |ψ t(e) E P 1 ∩(t(e)h(e)):···:P k ∩(t(e)h(e)) sq (N e (|ψ ψ| t(e) )), where if P j ∩ (t(e)h(e)) = / 0, we strip P j ∩ (t(e)h(e)) from the partition of E P 1 ∩(t(e)h(e)):···:P k ∩(t(e)h(e)) sq in the right-hand side, and |ψ t(e) is an arbitrary pure state which can be prepared at vertex t(e).
Capacities in this scenario are defined, similarly to Eqs. (3), (4), and (16)- (19). In particular, we just need to regard With this rephrasing, we can define various capacities for the multipartite communication.
As Eq. (15) gives upper bounds on the capacities in the scenarios for multiple-pair communication, Eq. (22) provides upper bounds on the capacities in the multipartite communication, through convex optimization (see Section II B 2). Lower bounds on the capacities can be obtained 69 by borrowing the idea of the aggregated quantum repeaters 63 , similarly to the scenarios for multiple-pair communication. In this case, by running entanglement generation over every quantum broadcast channel N e , the network shares a quantum state e∈Eρ e ε-close to e∈E (|GHZ 2 GHZ 2 | e ) ⊗l e R e δ associated with an undirected multi-hypergraph G = (V, E ), where each undi-rected hyperedge e ∈ E with the same ends of e corresponds to a GHZ state |GHZ 2 e among ends of e . Then, if we find a Steiner tree spanning clients C i , which is an acyclic subhypergraph connecting all vertices of C i , in the undirected multi-hypergraph G , we can provide a state close to a qubit GHZ state |GHZ 2 C i by performing generalized entanglement swapping 73 at vertices composing the Steiner tree. Hence, depending on an objective function in a definition of a capacity, by finding out edge-disjoint Steiner trees spanning clients C i in the undirected mutli-hypergraph G , we can give a state close toτ C d : However, this problem of finding the number of edge-disjoint Steiner trees, referred to as Steiner tree packing, is known to be an NP complete problem, even in the case of graphs 74 . Focusing on the scenario to distribute a single GHZ state or a single multipartite private state over a usual quantum network, Bäuml et al. provide 64 a computationally efficient method to obtain upper and lower bounds, which are explained in Section II D.
Recently, an alternative bound on the conference key and GHZ rates is provided by Das et al. 75 . In this work the entire network described by a so-called quantum multiplex channel, i.e. a quantum channel with an arbitrary number of in-and output parties and the protocol consists of many uses of the quantum multiplex channels interleaved by LOCC.

C. Flow problems in optimisation theory
A number of results provide upper and lower bounds on the rates at which entangled resources for the tasks described in Section II B can be distributed in a quantum network given by an arbitrary graph 32,42,43,45,63,64,68,69 . Besides, the bounds are closely related to standard problems in graph theory. Here we provide some background information on the graph and network theoretic tools, before seeing such a relation in detail (which are reviewed in Section II D). In particular, we review concepts such as network flows, the max-flow min-cut theorem, multicommodity flows, Steiner trees, as well as the complexity of the various optimisation problems occurring in this context. The discussion in this section applies to general flow networks in an abstract sense. Applications to quantum networks will follow in Section II D.
In the following we consider an undirected, weighted graph G = (V, E ), i.e. a graph consisting of vertices v ∈ V , undirected edges {vw} ∈ E where each edge is assigned a weight c {vw} ≥ 0. Note that in graph theoretic literature the weights are also referred to as capacities of the edges. Here we use the term weight when talking about abstract graphs to avoid confusion with information theoretic capacities in the following section.

Single commodity flows
We begin with the simplest scenario of a network flow, in which a single abstract commodity is assumed to flow through a network. For every edge {vw} ∈ E we can define edge-flows f vw ≥ 0 and f wv ≥ 0 such that The quantities f vw and f wv could be interpreted as the amount of the abstract commodity flowing from vertex v to vertex w and from vertex w to vertex v, respectively, while being constraint by a capacity c {vw} . Note that an alternative definition of flows which is commonly used is based on directed edges each of which is assigned a weight. Whereas the results we are going to present can be formulated in both scenarios, in our case it is convenient to use the version based on undirected graphs. Further, we choose two vertices s,t ∈ V , which we call the source and target, respectively. We can define a flow from s to t as while requiring that for every w ∈ V such that w = s and w = t it holds Eq. (24) can be regarded as the net flow leaving the source, whereas Eq. (25) can be regarded as flow conservation in every vertex except the source and target. By the flow conservation, it holds justifying the interpretation as a flow from s to t. The obvious task now is to maximise the flow f s→t defined by Eq. (24) with respect to the capacity constraint given by Eq. (23) and the flow conservation constraint given by Eq. (25), which is a linear program (LP). Namely ∀{vw} where the maximisation is over 2|E | edge flows f vw ≥ 0 and f wv ≥ 0. The LP given by Eq. (27) 77 .
Next, let us introduce the concept of cuts. Given a subset W ⊂ V of vertices, we define a W -cut as the set of edges i.e. as a set of edges the removal of which disconnects W and V \W . In the case where s ∈ W and t ∈ V \W , we call Eq. (28) an st-cut and use the notation ∂ (W s:t ). Further, we define the weight of a cut as For a given source and target, we can now define the minimum st-cut as We are now ready to introduce the max-flow min-cut theorem 78,79 , which states that for all s,t ∈ V it holds Let us note that so far we have only made the assumption that all edge flows are nonnegative real numbers. Depending on the network scenario considered, however, it might be necessary to require the values f vw and f wv to be integers for all {vw} ∈ E . In that case, the weights c {vw} can, without loss of generality, be taken to be integers as well. A graph with integer weights can equivalently be written as a multigraph, where each edge {vw} is replaced by c {vw} parallel edges, each with weight equal to one. The task of finding the maximum integer flow from s to t is then equivalent to finding the number of edge disjoint paths from s to t. It has been shown by Menger 66 , that the number of edge disjoint paths is equal to the minimum st-cut. Finally, let us note that adding an integer constraint turns the optimisation problem given by Eq. (27) into an integer linear program. Integer Linear programming has been shown to be an NP-complete problem 80 .

Multi commodity flows
The theory of network flows can be generalised to scenarios with k source-target pairs (s 1 ,t 1 ), ..., (s k ,t k ), and k different commodities flowing through a network concurrently. To describe such a scenario, for every edge, we define 2k edge flows f We can now define the flow of commodity i from s i to t i as Further we require that every commodity i is conserved in every vertex except its source s i and target t i , i.e. for all w ∈ V such that w = s i and w = t i it holds When maximising the flows we now have a number of possibilities for figures of merit. Here we focus on two commonly used figures of merit. The first option is to maximise the sum over the k flows, resulting in the maximum total concurrent flow where the maximisation is under the constraints that Eq. (32) holds for all edges {vw} ∈ E and Eq. (34) holds for all commodities i ∈ {1, ..., k} and for all vertices w ∈ V such that w = s i and w = t i . Eq. (35) is a linear program, which in a standard form has N = (2k + 1)|E | nonnegative variables and M = |E | + k(|V | − 2) equality constraints, and hence scales polynomially in the network parameters.
Whereas there might be situations where the total concurrent flow is a suitable figure of merit, we note that, depending on the network structure and weights, Eq. (35) could be maximised by a flow instance where the f s i →t i vary greatly. The second figure of merit we consider here is fairer in the sense that it maximises the least flow which can be achieved for any of the commodities concurrently. Namely, we have the worst case concurrent flow where again the maximisation is under the constraints that Eq. (32) holds for all edges {vw} ∈ E and Eq. (34) holds for all commodities i ∈ {1, ..., k} and for all vertices w ∈ V such that w = s i and w = t i . Adding a slack variable f worst ≥ 0, Eq. (36) can be reformulated as the maximisation of f worst such that, in addition to the above constraints, for all i ∈ {1, ..., k} it holds f worst ≤ f s i →t i , which is a linear program. This in a standard form has N = (2k + 1)|E | + 1 + k nonnegative variables and M = |E | + k(|V | − 2) + k equality constraints, again scaling polynomially in the network parameters. Let us note that both multicommodity flow maximisations, given by Eqs. (35) and (36) reduce to the flow maximisation given by Eq. (27) in the case of a single commodity, i.e. k = 1. An obvious question arising now is whether the maxflow min-cut theorem, Eq. (31) can be extended to the case of multicommodity flows. Again, there are different quantities to consider. The first way to generalise the minimum st-cut to multiple source-target pairs is to consider multicuts. For k commodities, (s 1 : t 1 ), ..., (s k : t k )-multicut is a set of edges the removal of which disconnects s i from t i for all i ∈ {1, ..., k}. As in the case of cuts, the weight of a multicut is defined as the sum of weights of the edges it contains. Minimising over all multicuts, we obtain the minimum multicut c multicut min . In the case where k = 1, the problem of computing the minimum multicut reduces to the minimum cut, and hence it can be solved in polynomial time via the max-flow min-cut theorem. Surprisingly, even for k = 2, the problem can be reduced to a single-commodity flow 81 and can be computed in polynomial time. However, in general, the problem of finding the minimum multicut is known to be NP-hard 82 .
A second way in which the minimum st-cut can be generalised to multiple source-target pairs is the min-cut ratio, which is defined as where the minimisation is over arbitrary subsets W ⊂ V and d(∂ (W )) is the number of source-target pairs separated by ∂W . Computation of the min-cut ratio has been shown to be NP-hard in general, as well 83 .
Whereas it is clear from a complexity point-of-view that, for general k, there cannot be an equality relation between c multicut min or c cut ratio min to a multi-commodity flow instance that can be computed in polynomial time, there are relations up to a factor of O(log k). Namely, it has been shown by Garg et al. 84 that where Similarly, it has shown by Aumann et al. 83 that where

Steiner trees
When considering scenarios involving groups consisting of more than two network users rather than pairs of users, we need to generalise the concept of a flow from a source to a target or, in the discrete case, the concept of edge-disjoint paths in a unit-capacity multigraph to a multipartite setting. In the discrete case, this leads to the concept of Steiner trees. Namely, given a unit-capacity multigraph G = (V, E ) and a set S ⊂ V of vertices, a Steiner tree spanning S, or in short S-tree, is a connected, acyclic subgraph of G connecting all s i ∈ S. By definition, in the S-tree, any pair of vertices s i , s j ∈ S with i = j is connected by exactly one path. In the special case where S = V , an S-tree is also a spanning tree. The relevant task now is to find the number of edge disjoint S-trees in G . This problem, which is known as Steiner tree packing, is NP-complete in general 85 . Note however that in the case where |S| = 2, the problem reduces to finding the number of edge-disjoint paths, which by Menger's theorem is equal to the min-cut.
In the general case, it is possible to upper bound the number of edge-disjoint S-trees by the S-connectivity of the graph G : As in an S-tree any pair s i , s j ∈ S with i = j of vertices has to be connected by a path, the number t S (G ) of edge-disjoint S-trees in G cannot exceed the number of edge-disjoint paths between any pair s i , s j ∈ S with i = j in G . Minimised over all s i , s j ∈ S with i = j, this is known as the S-connectivity λ S (G ) of the graph. Application of Menger's theorem to each pair s i , s j ∈ S with i = j individually shows that λ S (G ) is equal to the minimum S-cut, which is defined as Further, there have been a number of results [86][87][88] providing lower bounds on the number of edge disjoint S-trees t S (G ), which are of the form Kriesell 86 conjectured that g 3 = 1 2 and g 4 = 0. Lau 87 has shown that Eq. (41) holds for g 3 = 1 26 and g 4 = 0. Finally, Petingi et al. 88 show the relation for g 3 = |V \S| 2 and g 4 = 1.

D. Efficient bounds for communication tasks
As was shown by Bäuml et al. 64 , a combination of the results presented in the two previous sections provides efficiently computable upper and lower bounds on the rates at which various resource states can be distributed in a network. The upper and lower bounds are in the form of single-and multi-commodity flow optimisations in an undirected version of the original network graph. The weights of the corresponding edges are in terms of the quantum capacities Q (for lower bounds) and the entanglement E of the channels (for upper bounds). Thus, given one knows these quantities for every channel constituting the network, one can obtain bounds on the network capacities by simply running linear programs.
Given a quantum network described by a directed graph G = (V, E) as defined in Section II B, we define an undirected graph G = (V, E ) as follows: if there is only a single directed edge connecting two vertices v, w ∈ V , we regard the edge as an undirected edge of the graph G with the same two ends; if there exist two edges vw ∈ E and wv ∈ E for two vertices v, w ∈ V , we only add one undirected edge {v, w} = {w, v} to E . The reason we switch from the directed graphs to undirected ones here is that, whereas the original quantum channels N e for e ∈ E are directed, in protocols such as the aggregated repeater one presented in Section II B, those channels are used to distribute Bell pairs, which are intrinsically undirected under unlimited LOCC, in the beginning. The remainder of the protocol is then concerned with transforming the resulting Bell-state network into the desired target state.
Next, each edge {vw} ∈ E of G can be equipped with two weights c E ({vw}) and c Q ({vw}) to be used in the flow maximisations for the upper and lower bounds, respectively. For protocols in Λ(n, {q e } e∈E , ε), we define where we set E (N vw ) = Q(N vw ) = 0 whenever vw / ∈ E.

Bipartite communication
We begin by considering the bipartite case where two vertices, A and B, wish to obtain Bell states or bipartite private states. The corresponding capacities are given by Eqs. (3) and (4). Looking at Eq. (5), we now see that the upper bound presented there is indeed equal to the min-cut in the undirected graph G with weights c E ({vw}) for all {vw} ∈ E . By the max-flow min-cut theorem, Eq. (31), it can be expressed as a flow maximisation of the form Eq. (27), which can be computed by linear programming. As the constraints in the maximisation over usage frequencies q e in Eq. (6) are linear as well, we can now formulate the entire upper bound given by Eq. (6) as a linear program.
Similarly, the lower bound on the capacity provided by Eq. (14) can be transformed into a flow maximisation in the undirected graph G with weights c Q ({vw}), and hence the optimisation in Eq. (14) becomes a linear program. In summary, we have 64 max where we have defined f A→B Q max (G ) and f A→B E max (G ) as maximum flow from A to B in G with edge weights given by Eq. (43) and (42), respectively. Further, the maximisation is over q e ≥ 0 such that ∑ e∈E q e = 1 on both sides. Both the upper and lower bounds are linear programs of the form (27), with additional optimisation over usage frequencies, which adds |E| non-negative variables and one equality constraint, and hence still computable in polynomial time.
The reason we can drop the integer constraint in the lower bound is that we are looking at the asymptotic limit 64 . Roughly speaking, assuming that the f A→B Q max (G ) is maximised by non-integer edge flows f vw , we can always find a large enough number N such that N f vw are close enough to an integer for all edges. In particular, by using the corresponding channel a large enough number of times, we can obtain N f vw Bell pairs for all edges {vw} ∈ E , which can then be connected by entanglement swapping along paths linking Alice and Bob.

Multi-pair communication
Let us consider the case of multiple user pairs (A 1 , B 1 ), ..., (A k , B k ) wishing to obtain bipartite resources such as Bell pairs or private states between them concurrently. In Section II B it was shown that there are a number of different ways to define the corresponding capacities. One is to maximise the sum of the rates achievable by the user pairs. The corresponding quantum and private multi-pair network capacities Q total and P total , respectively, are defined in Eq. (19) with s i = 1 for all i. A second way is to maximise the worst case rate achievable between all user pairs, with corresponding network capacities Q worst and P worst , defined in Eq. (17).
It follows as a special case from the results of Bäuml and Azuma 69 that the total private multi-pair capacity P total is upper bounded by the minimum multicut c multicut E sq min (G ). Similarly, it has been shown 64,69 that the worst-case private multi-pair capacity P total is upper bounded by the minimum cut ratio c cut ratio E min (G ) with respect to weights defined by Eq. (42) for all entanglement measures E holding the bound (1). Applying the generalisations of the max-flow min-cut theorem to multi-commodity flows, given by Eqs. (38) and (39), respectively, we can upper bound both private multi-pair capacities in terms of multi-commodity flow optimisations. Using a generalisation of the aggregated repeater protocol 63 , lower bounds on the multi-pair quantum capacities Q total and Q worst can be achieved by the corresponding multi-commodity flow maximisations f total Q max (G ) and f worst Q max (G ) with respect to weights defined by Eq. (43) 64 . Again we note that the integer constraint in the flow optimisation can be dropped in the asymptotic case. In summary, it holds for the total multi-pair network capacities where the maximisation is over all q e ≥ 0 with ∑ e q e = 1. Both upper and lower bounds are linear programs. Similarly, for the worst-case multi-pair network capacities can be upper and lower bounded by polynomial-size linear programs, as where E can be any entanglement measure for which the bound (1) holds. Again, the upper and lower bounds are polynomial-size linear programs. As we discuss in the next section, the gaps g 1 (k) and g 2 (k), both of which are of order O(log k), leave the possibility for more efficient protocols, e.g., quantum network coding. Let us also note that the multi-commodity flow based approach to obtain achievable rates in multi-pair communication can be extended to realistic scenarios with imperfect swapping operations and without entanglement distillation. In such a case, the fidelity achievable between the user pairs drops with the number of swapping operations. Chakraborty et al. 89 overcome this problem by adding an additional constraint on the lengths of paths to the optimisation of the total multicommodity flow. They provide a linear program scaling polynomially with the graph parameters, which provides a maximal solution if it exists.

Multipartite communication
Finally we consider the case where a set of users A = {A 1 , ..., A k } wishes to establish a k-partite GHZ or private state. It is straightforward to generalise the bipartite quantum and private network capacities given by Eq. (4) to the multipartite case. We denote the corresponding multipartite quantum and private network capacities Q A (G, {N e } e∈E ) and P A (G, {N e } e∈E ), respectively. From the results of Bäuml and Azuma 69 it follows that the multipartite private network capacity is upper bounded by the minimum A -cut in G with weights c E sq ({vw}): By application of the max-flow min-cut theorem (for every pair A i , A j ∈ A with i = j), as well as introduction of a slack variable, the r.h.s. can be expressed as a linear program 64 .
Namely we maximise f ≥ 0 such that for all pairs A i , A j ∈ A , i = j, f ≤ f A i →A j and requiring that the capacity and flow conservation constraints, given by Eq. (23) and (25), respectively, are fulfilled for every pair.
In order to lower bound Q A (G, {N e } e∈E ), we have to generalise the aggregated repeater protocol. The first part of the protocol, in which Bell states are distributed resulting in a Bell network described by a multigraph G , remains the same. The second part, however, has to be generalised to the distribution of GHZ states among the parties in A . Namely, instead of connecting paths of Bell pairs by entanglement swapping, it is possible to connect A -trees of Bell pairs by means of a generalised entanglement swapping protocol 73  Whereas the number of edge disjoint A -trees is an NP-hard problem, we can make use of Eq. (41) and Menger's theorem to relax the problem an integer-flow optimisation, which in the asymptotic limit can be further relaxed to a flow optimisation problem of the form min A i ,A j ∈A ,i = j f A i →A j Q max , which is a linear program. In summary we have the following bounds: where we have set g 3 = 1 2 . Both the upper and lower bounds are linear programs with 2 |A | 2 |E | + |E| + 1 variables and |A | 2 (|E | + 1) + |E| inequality constraints, i.e., of polynomial size.

E. Beyond classical routing in quantum communication networks
So far we have focused on the distribution of Bell, GHZ and private states by means of protocols based on the aggregated repeater protocol, which basically is a classical routing protocol applied on quantum networks. Other quantum network protocols based on classical routing have been introduced, among others, by Pant et al. 91 , Schoute et al. 92 and Chakraborty et al. 93 . We will now give a brief, and by no means exhaustive, overview of other quantum network protocols, which go beyond the classical routing paradigm.

Quantum network coding
The lower bounds presented in the previous section are based on routing protocols, where Bell pairs are established between nodes connected by channels and then connected by means of entanglement swapping along paths connecting the users 63 . Whereas, for a large class of quantum channels, in the single pair case, the upper bounds can be achieved in this way, the situation becomes more involved in the multipair case, where we have seen gaps between upper and lower bounds. Such gaps leave room for improvement in terms of more efficient protocols such as a quantum version of network coding 94 .
In classical network theory, network coding protocols differ from routing protocols in that bits of information are being acted upon by some operation before being sent via a channel, allowing the combination of several bits of information into one. Whereas routing is optimal in single-pair scenarios, coding can provide an advantage in multi-pair problems in directed networks, as can be seen in the famous butterfly graph. In undirected networks, the question of a coding advantage in the multiple unicast setting is still an open question 95 . In the case of directed classical multicast settings, where messages are sent from a source node to several destination nodes, it has been shown that coding can always achieve the capacity, which is equal to the so-called cutset bound, whereas routing cannot, which can, again, be demonstrated using a butterfly graph 96 .
Returning to the setting of quantum networks, a quantum version of classical network coding has been shown to outperform the simple routing approach presented above 97,98 . Namely, assuming we have a butterfly network where each link consists of a single Bell pair, concurrent quantum routing between the two diagonal node pairs is not possible as any path connecting the first pair of diagonal nodes disconnects the second pair when being used up. However, using the Bell states as identity channels via teleportation, we can apply the network coding protocol presented by Kobayashi et al. 97,98 , to distribute Bell pairs concurrently between both user pairs. The protocol mainly consists of creation of the state |+ = (|0 +|1 )/ √ 2 at each node followed by a translation of a classical network coding protocol into sequences of Clifford operations. The resulting encoded states are transmitted via the identity channels. A similar advantage in a butterfly network over the Steiner tree routing approach can be obtained when distributing GHZ states using quantum routing 99 .
Going to the asymptotic setting, on the other hand, routing in combination with time-sharing, i.e. the distribution of entanglement between one pair in one round and another pair in another round and so on, is sufficient to achieve the cut based upper bounds in the multi-pair scenario 100 . Note that there is a tradeoff between the increased need of quantum memory in time-sharing and the need to perform more complex operations in quantum network coding.

Graph state protocols
The aggregated repeater protocol discussed above consists of an initial stage distributing Bell states across all edges of the graph, resulting in a Bell state network involving all connected vertices in the graph. This Bell state network is then transformed into the desired target state among a single or multiple pairs or groups of users. Thus, the Bell state network can be seen as a universal resource state for a number of network tasks. A different strategy is to create a large multipartite entangled graph state among all, or a large subset of, the vertices, which serves as a universal resource state that can subsequently be transformed into the desired target states. In this class of protocols there are two tasks to be considered. The first is to use a network of quantum channels and free operations (such as LOCC) to create the graph state that serves as resource. The second is to transform the graph state into the desired target states by means of free operations.
The first task has been considered by Epping et al. 101 , where the goal is a multipartite entangled graph state among all vertices in the network, except the ones that merely serve as repeater stations to bridge long distances. This is achieved by protocol involving local creation of qubits in |+ -states, controlled-Z operations entangling two qubits, Pauli-X measurements and the sending of qubits via channels, corresponding to directed edges. In a first step all vertices are connected to form a graph state, which is done as follows: In vertices with only n outgoing edges, n + 1 |+ -states are created, all of which get entangled by means of controlled-Z operations and n qubits are sent via the outgoing edges. In vertices with incoming and outgoing edges, all incoming qubits get entangled among each other as well as with locally created qubits in |+ -states which are then sent via the outgoing edges. In a second step the repeater stations are being removed from the graph state by Pauli-X measurements. This protocol is then combined with stabiliser codes to compensate for losses in the transmission channels, noisy gates as well as errors in the preparation and measurement. In subsequent work 102 the technique is combined with quantum network coding and generalised to qudit graph states.
The second task has been considered by Hahn et al. 103 , which is concerned with the transformation of a graph state into Bell states between a single or multiple pairs of users, GHZ states among groups of users and, most generally, the transformation from one graph state to another graph state. The free operations are restricted to local Clifford operations and Pauli measurements. Whereas it is always possible to isolate the shortest path connecting two users and make the connection by measuring the intermediate vertices, Hahn et al. provide a protocol that requires the measurement of less nodes. Their approach is based on a graph transformation known as local complementation 104 and the fact that a graph state can be transformed into another graph state by local Clifford gates whenever the corresponding graphs can be interconverted by means of local complementations 105 . In related work by Dahlberg et al. 106,107 , the question of interconvertability between graph states by means of single-qubit Clifford operations, single-qubit Pauli measurements and classical communication is studied. The authors show that the problem of deciding whether two graph states can be interconverted by this class of free operations is NP-complete and present an algorithm that provides the sequence of operations necessary to perform the task, if possible. Another approach is presented by Pirker et al. 17 , where tensor products of GHZ states or fully connected decorated graph states are used as resource states, allowing for creation of arbitrary graph states among the network users.

Coherent control
The adaptive protocols described in Section II B, while using quantum channels, quantum states, and quantum de-and encoding operations, are classical in the way they choose which channel is used in each round: One could think of the LOCC operations having a classical output register that is used as a control register determining which channel is used next. There exist, however, scenarios in which the choice of channels is not determined by a classical register, but coherently using a quantum control register. For example, a 'quantum switch' has been introduced, where a quantum control register determines the order of two channels N 1 and N 2 . By entering a state in quantum superposition, say |+ c into the control register, a superposition of the orders of N 1 and N 2 is achieved 108 .
Surprisingly, it has been shown that such a superposition of orders can boost the rate of classical and quantum communication beyond the limits of conventional quantum information theory [109][110][111] . In particular, it has been shown that two identical copies of a completely depolarising channel can become able to transmit classical information when used in such a superposition of orders. Similar effects were also observed when quantum control is used to create a superposition between single uses of two channels 112 .
Miguel-Ramiro et al. 113 have introduced the concept of coherent quantum control into the context of quantum networks, allowing for coherent superpositions of network tasks. In particular they provide protocols where two nodes in a network are connected by a superposition of different paths. Other examples include the distribution of a quantum state to different destinations in superpositions and the distribution of a superposition of GHZ states among different user groups. We note that such protocols are not included in the paradigm of adaptive LOCC protocols used in Section II B.

III. WAITING TIMES AND FIDELITY ESTIMATION OVER ABSTRACT QUANTUM NETWORKS
In this section, we review analytical tools for characterizing the performance of quantum networks and the algorithms that immediately follow the analytical expressions. In particular, we consider the literature that studies the time it takes to distribute remote entanglement over a quantum network, referred to as the waiting time, and the quality of the entanglement. Due to their more modest quantum information processing requirements, we devote a large part of the section to quantum repeaters which are built from probabilistic schemes, i.e. the so-called first-generation repeater 114 . As a consequence of the probabilistic nature of such schemes, the waiting time is a random variable; thus, it is not represented by a single number but instead by a probability distribution. Our presentation focuses on the fidelity with respect to the desired maximally entangled state as a measure of entanglement quality. However, many of the tools can directly be used for estimating other figures of merit such as the secret key rate.
This section is organized as follows. We start in Section III A with the modeling of a quantum network, which includes the mathematical abstraction of different components in a quantum network. In Section III B, we discuss the analytical tools used in evaluating the performance of networks. In some cases, those analytical tools yield closed-form expressions, of which the evaluation requires the assistance of numerical algorithms. We discuss three such cases in Section III C: Markov chain methods, probability-tracking algorithms, and sampling with Monte Carlo methods. Finally, we review in Section III D the analysis of quantum repeater protocols that include quantum error correction. We have chosen to limit the scope of this section to discrete variable protocols.

A. Abstract models of quantum networks
Here, we summarize common models of different network components, with an emphasis on how they contribute to the statistics of the waiting time and quality of the entangled state. Throughout the section, we refer to a pair of entangled qubits shared by spatially-separated nodes, as a 'link'.
Entanglement generation. By entanglement generation, we refer to the delivery of a fresh Bell state between two nodes in the network which are directly connected through a communication channel, such as an optical fiber. We refer to the generated entanglement as an elementary link. There are several schemes for the generation of elementary links 114 , and in each of them, the generation is performed in discrete attempts until the first successful attempt. We assume that each attempt is of constant duration ∆ and has constant success probability p g . The attempt duration ∆ is determined by the distance and speed of light in the medium; in the rest of this section, we set ∆ = 1 for simplicity. It is also commonly assumed that the distinct attempts are independent and thus that the state ρ that is produced is constant, i.e. it is independent of the number of attempts required to produce it. The state ρ is a noisy Bell state which typically incorporates different sources of noise, photon loss, and detector inefficiency.
Entanglement swapping. The photon transmission probability in fiber decays exponentially and thus fundamentally limits that distance over which elementary link generation can be performed 44 . This limit can be overcome by adding additional nodes in between, so-called quantum repeaters, which perform entanglement swaps to connect two short-distance links into a single long-distance one 115 . Typically, entanglement swaps are probabilistic, with a fixed success probability p s which is normally independent of the states swapped but depends on the physical implementation. For matter qubits that can be controlled directly, an entanglement swap can be implemented with deterministic quantum gates, i.e. p s = 1. If entanglement swapping is implemented with optical components, the entanglement swapping becomes probabilistic, i.e., p s < 1 and typically p s ≤ 0.5 116 . There are also more sophisti-cated optical swapping schemes with a probability larger than one half [117][118][119] . In some models, where the memory decoherence to the vacuum state is considered, the success probability can also be a variable 120 .
Entanglement distillation. Performing an entanglement swap on two imperfect Bell states yields long-distance entanglement of lower quality than the original two states individually had, and its quality is further decreased by imperfections in the quantum operations. In principle, entanglement distillation can be used to improve the fidelity by probabilistically converting multiple low-quality entangled pairs of qubits into a single one of high quality 121 . In contrast with entanglement swapping, the success probability p d depends on the entangled states that are distilled 121,122 .
Entangled state representation. Arguably, the simplest model of the fresh elementary link state is a Werner state 123 , which characterizes the state with a single parameter w: where |Φ 2 is the desired maximally-entangled two-qubit state andÎ/4 the maximally mixed state on two qubits. Although operations such as entanglement distillation do not always output a Werner state, any two-qubit state can be transformed into a Werner state with LOCC without changing the fidelity 25 . A more general model is a probabilistic mixture of the four Bell states. This representation is convenient as it includes the resulting state after the application of random Pauli gates on a perfect Bell state. In principle, one could also track the full density matrix, though many studies choose the previous two representations to simplify the analysis. Given the density matrixρ of a state, its fidelity with a pure target state |φ is given by φ |ρ|φ . Throughout the section, the target state will be a Bell state.
Noise modeling. Imperfections of the quantum devices, for example, operational noise and detector inefficiencies, are commonly modeled by depolarizing, dephasing, or amplitude damping channels. The first two can be incorporated relatively simply into analytical derivations as they correspond to the random application of Pauli gates. Amplitude damping requires tracking the full density matrix. One could, however, replace an amplitude damping channel with the more pessimistic choice of a depolarization channel, which does not change the output state's fidelity with the target state, or alternatively twirl the damped state by applying random Pauli operations 124 .
Particularly relevant in the context of entanglement generation using probabilistic components is the noise caused by time-dependent memory decoherence: in case multiple links are needed, the earliest link is generally generated before the others are ready and thus needs to be stored in a quantum memory. The storage time leads to a decrease in the quality of the entanglement, and the longer the qubit is stored, the more its quality degrades. Due to the interplay between waiting time and time-dependent decay of entanglement quality, memory noise is particularly hard to capture. Sometimes this problem is sidestepped by analyzing protocols with running time qualitatively shorter than the memory decoherence time.
Node model. For simplicity, the network nodes can be modeled by a fully-connected quantum information processing device capable of generating entanglement in parallel with its neighbors. However, it is important to note that many platforms do not conform to this model. For instance, NV-centers in a single diamond have a single optical interface. Hence, if nodes hold only a single NV center, entanglement generation can only be attempted with one adjacent node at a time. Moreover, the connectivity between the qubits follows a star topology, i.e. direct two-qubit gates between arbitrary qubits are not possible.
Cut-off. Due to memory decoherence, the quality of the stored entanglement decreases as the waiting time grows. One common strategy to compensate for memory decoherence is cut-offs: if a link remains idle for too long, it is discarded. By discarding entanglement whose storage time exceeds some pre-specified threshold, one improves the quality of the delivered entanglement at the cost of longer waiting time.
Additionally, it is possible to build on top of this idea a simplified model of memory decoherence: the quantum information is preserved perfectly for a fixed duration and then lost 125,126 .
Note that the inclusion of cut-offs in entanglement distribution schemes complicates their analysis because of the additional effect of waiting time on the state quality.
Nested protocols. One particularly relevant network topology is the repeater chain, where all nodes are arranged in a line. Nested protocols offer a structured approach to distributing entanglement across a repeater chain [127][128][129][130][131][132][133][134] . In this section, unless explicitly mentioned, we follow the BDCZ scheme 127 , with the restriction that each entanglement swap doubles the distance that an entangled pair spans. In such a scheme, the number of repeater segments is 2 n (2 n + 1 nodes) where n is the number of nesting levels at which an entanglement swap is performed. We depict examples of BDCZ protocols in Fig. 2. We denote the waiting time random variable of a repeater scheme on 2 n segments as T n . Many of the tools for determining the waiting time statistics and quality of the produced entanglement discussed below also apply to other schemes than nested repeater protocols.

B. Analytical study of the waiting time and fidelity
In this section, we present analytical tools to compute the waiting time and the fidelity of the entangled state produced between the end nodes of a repeater chain.
We consider the nested repeater chain protocol on 2 n segments (see Section III A) with only entanglement swapping, i.e. no distillation or cut-offs unless explicitly mentioned. For simplicity, we assume that the generation probability p g is the same for each pair of adjacent nodes and the swapping probability p s is equal at each nesting level. We also assume that the nodes are capable of generating entanglement in parallel. Finally, we ignore the (constant) duration of local operations and classical communication for simplicity, although all of the tools mentioned are capable of incorporating these.
We first investigate methods that instead of tracking the full At each nesting level, the distance that entanglement spans is doubled. By T n , we denote the random variable describing the delivery time of entanglement at level n. In Section III, we consider nested repeater protocols of entanglement generation and swapping such as in (b). Whenever the protocol includes entanglement distillation, as in (a), or cut-offs, it is mentioned explicitly.
probability distribution, only track an approximation of the average waiting time and quantum state at each nesting level of the entanglement-distribution protocol. To demonstrate the tools used in computing the distributions, we include an explicit calculation for a protocol on two nodes with a single repeater positioned in between. As this exact calculation cannot be directly generalized to a higher nesting level in a nested protocol with more than a single repeater, we then consider the idea of approximating the waiting time by assuming it follows the statistics of elementary-link generation, where the mean waiting time is computed using the approximation from the previous level. Finally, we review the mathematical tools and approximation methods used to analyze deterministic swapping protocols and distillation-based repeater schemes.

The mean-only approach
In many early analyses of repeater protocols, only the mean waiting time is considered for each nesting level: it is assumed that the entanglement is delivered after a fixed time duration determined by the generation rate. We refer to it as the mean-only approach. In this approach, the mean waiting time is computed as the product of the mean waiting time at each nesting level (1/p g at the bottom level, 1/p s at the higher levels), yielding T n = 1/p n s p g 127, 128 . This approach can be refined by noting that at each nesting level the protocol proceeds only when two adjacent pairs are ready. Then, the mean waiting time can be approximated by T n = (3/2) n /p n s p g 130,133,135-140 . The factor 3/2 comes from the fact that in the limit of very small success probability, the waiting time of preparing two links is approximately 3/2 times that of one link. We discuss the statistics behind this factor later in Section III B 3.
The mean-only approach is a good approximation when p g and p s are much smaller than 1 135,136 . However, it only approximates the mean, i.e. it does not provide the entire probability distribution of the waiting time. Hence, it is not suited for investigating time-dependent aspects such as memory noise or cut-offs. With this method, memory decoherence is either approximated by an inefficiency constant 139 or studied only for the communication time 141 . To provide a better estimation, one needs to consider the waiting time distribution and the statistics it results in, which we discuss below.

Single repeater swap protocol
Here, we explicitly compute the probability distribution of the waiting time for the simplest repeater chain: a single repeater between two end nodes. We also derive an expression for the mean fidelity decay due to memory decoherence. Many problems regarding single-repeater protocols have an analytical solution because the entanglement generation follows a known distribution. By studying this simple scenario, we demonstrate the common concepts and methods used to describe and compute the statistics of waiting time and fidelity. In later sections, we use this calculation as a basis for the analysis of nested repeater protocols of more nodes.
We describe the waiting time of elementary-link generation as a random variable T 0 , following a geometric distribution given by where t ∈ {1, 2, 3...} and p = p g . This distribution plays a central role in the statistics of entanglement distribution, as we see in the remaining part of this section. In the limit of p g 1, the geometric distribution can be approximated by an exponential distribution, which is a continuous distribution with t ≥ 0. Note that we have set the attempt duration ∆ of entanglement generation to 1 (Section III A).
To perform an entanglement swap, both elementary links have to be prepared first. We define the time used in preparing them as M 0 : where T 0 is an independent copy of T 0 . The distribution of M 0 can be computed using the fact that for any independent random variables X and Y . The mean of M 0 , i.e. the waiting time until both elementary links have been prepared, is given by 137 : After two elementary links are prepared, the repeater node performs an entanglement swap, which is a probabilistic operation with success probability p s . The total waiting time for generating the entanglement between two end nodes is therefore where K represents the number of swap attempts until it succeeds and M (k) 0 are independent copies of M 0 . Eq. (54) is referred to as a compound distribution since the number of summands K is also a random variable. For an entanglement swap, the number of attempts K also follows a geometric distribution (Eq. (49)) with success probability p = p s . Because K and M 0 are independent, the average waiting time is given by The intuition behind Eq. (55) is that, on average, the repeater node requires K swap attempts until the first successful swap, and for each swap attempt, M 0 attempts are needed to prepare the two elementary links. Computing the fidelity of the two elementary links just before swapping can be done as follows. If the generation of elementary links is not deterministic, i.e. if p g < 1, the two elementary links are in general not produced at the same time, requiring the earlier of the two to be stored in a quantum memory. This storage time results in decoherence of the earlier link. To estimate the fidelity decrease, we need to first compute the distribution of the memory storage time, i.e. the time difference between the generation of the earlier and the later link. We define q g = 1 − p g . The probability that one link is prepared j steps before the other is given by 138,142 Pr Here we assume that T 0 > T 0 . Modelling the fidelity decrease as exponential-decaying function of the storage time, the fidelity of the earlier link decays by a factor Γ = exp(−|T 0 − T 0 |/t coh ) , where t coh denotes the memory coherence time. Plugging in Eq. (56), we obtain The factor 2 before the sum corresponds to the possibility that either link can be generated earlier than the other. Finally, the evaluation of the sum gives 138,142 .
In addition to the single-repeater scenario considered above, analytical results for the memory decay have also been obtained for more advanced single repeater protocols such as a protocol with cut-offs 125 or protocols where two elementary links have to be prepared sequentially 143 .
Unfortunately, for higher-level nested protocols, i.e. n ≥ 1, there is no analytical expression known for the mean waiting time T n with p s < 1, because T i for i > 0 does not follow a geometric distribution, in contrast to T 0 .

Approximation with geometric distribution at higher level
Above, we computed the waiting time probability distribution in the single-repeater scenario. This calculation explicitly relied on the fact that the waiting time distribution of elementary-link generation follows a simple distribution, the geometric distribution (Eq. (49)). Unfortunately, for nested repeater chains with more than a single repeater, no exact expression for the waiting time distribution has been found.
However, the waiting time distribution at higher nesting levels can be approximated by assuming that, at a higher level, the waiting time distribution is still geometrically or exponentially distributed (Eq. (50)). This approximation is usually used in an iterative manner. One computes the average waiting time at the current level and uses it to define a geometric distribution with the same expectation value. This new distribution is then used to study the next nesting level. In Fig. 3, we compare the approximated distribution with the exact one.
Let us give the explicit calculation under the approximation that the waiting time follows a geometric distribution at each nesting level. We first calculate T n−1 and then approximate the distribution of T n−1 with a geometric distribution (Eq. (49)) parameterized by p = 1/ T n−1 . Under this assumption, the mean waiting time T n can be computed by  3. The probability distribution of the exact waiting time T 2 of a nested swap protocol with 4 repeater segments (computed with the algorithm from Section III C 2) and, as an approximation to the exact distribution, the geometric distribution from Eq. (49) with the same mean, i.e. p = 1/ T 2 . (Top) We see that the two distributions deviate most for short waiting times. This can be explained by noting that the exact probability that all entanglement generation steps and entanglement swaps succeed in the first few steps is very small. This fact is not captured by the approximation. (Bottom) We observe that for small swap success probabilities p s (both axes are rescaled to compare only the shape of the distributions), the deviation becomes smaller.
a derivation analogous to the one leading to Eq. (55) in Section III B 2 and is given by with p n−1 = 1/ T n−1 . In the limit of p g → 0 for the bottom level (n = 0) and p s → 0 for higher levels (n > 0), the mean waiting time T n−1 → ∞ and thus p n−1 → 0. As a consequence, Eq. (58) can be approximated as Effectively, it means that, on average, the waiting time of generating two links is approximately 3/2 times that of a single link. Applying Eq. (59) iteratively over all nesting levels with T 0 = 1/p g yields which is precisely the 3-over-2 approximation mentioned in Section III B 1. The error introduced by the approximations Eq. (60) and Eq. (58) is shown in Fig. 4. As expected, the figure shows that the approximations behave well if the success probabilities of elementary-link generation and swapping are small, i.e. p g → 0 and p s → 0. The figure also shows that the approximations are not so good for large success probabilities; the deviation from the exact mean waiting time increases as p s grows, and the deviation is worse for Eq. (60) than for Eq. (58).
To approximate the fidelity of the produced link between the end nodes of the repeater chain in the presence of memory decoherence, one can use Eq. (57) by replacing p g with 1/ T n−1 . The approximation is computed analogously to the procedure described for the average waiting time; that is, for a given level, the average infidelity for the entanglement links just before a swap due to the memory decoherece is calculated and used to derive the initial infidelity for entangled links at the next level. By assuming that the distribution at every level is given by the exponential distribution (Eq. (50)), Kuzmin et al. designed a semi-analytical method for computing fidelity with more sophisticated memory decay models 145,146 (see also Section III C 2).
A different approach to keep the waiting time distribution geometric at a higher level is to design a special protocol. For example, Santra et al. 147 introduce a family of protocols with a memory buffer time. This buffer time is a threshold on the total waiting time for preparing the two links for the swap. If the links are not ready before the buffer time is reached, the protocol aborts and restarts from entanglement generation. The buffer time is slightly different from the memory cut-off (see Section III A); with a memory cut-off the protocol aborts if the memory storage time of a single link (instead of both links) exceeds a threshold.
The protocol is designed such that the buffer time at the current level becomes the time step at the next level. As a consequence, the waiting time is geometrically distributed at each nesting level. Note that the protocol results in avoidable additional memory decay as both links have to wait until the buffer time is reached even if they are prepared before that. Despite this, by optimizing the buffer time, Santra et al. show that this protocol improves the final fidelity for some parameter regimes compared to the nested repeater protocol without buffer times.
An alternative approach was taken for optimizing repeater protocols including distillation, where the buffer time is chosen large enough so that the protocol becomes nearly deterministic 148 . At a cost of longer waiting time and lower fidelity, the variance in the fidelity is kept small and the protocol can deliver entanglement at a pre-specified time with high probability.

Deterministic entanglement swap
So far we have focused on the regime where the success probability of entanglement swap is smaller than 1. In this section, we discuss nested repeater protocols with deterministic entanglement swapping (p s = 1) and without distillation.
First, let us compute the waiting time probability distribution in the case of deterministic swaps without distillation or cut-offs. Recall that we assume that the time required to perform local operations is negligible so that the deterministic entanglement swap has no contribution to the waiting time. We plot the relative difference between the approximated and the exact mean waiting time as a function of the swap success probability p s for nested swap protocols up to level n = 4. The three approximation methods are: 1. The approximation with geometric distribution given by Eq. (58). It is exact at level 1 and the deviation increases as the levels grow; 2. The 3-over-2 approximation given by Eq. (60), which is itself an approximation of Eq. (58). In the limit of p = p g 1, the difference between the two approximation vanishes; 3. The deterministic swap approximation (Eq. (63)) assuming that the swap always succeeds (p s = 1). Note that, in contrast to the former two, the latter approximation is a lower bound on the exact waiting time. The deviation of the deterministic swap approximation is almost independent of p g . The exact mean waiting time T E is computed using the algorithm in Section III C 2. Due to our implementation, we cannot reach the region p s → 0. However, Shchukin et al. numerically show that, in this limit, the deviation of the approximation with geometric distribution becomes negligible (described in Ref. 144 ).
i.e. the maximum of N independent copies of T 0 : The cumulative distribution of T N from Eq. (61) is given by the general version of Eq. (52) for the maximum of N independent and identically distributed random variables: from which the probability distribution of T N can be computed using By T N,k , we denote the random variable describing the time at which the first k elementary links of the N segments are generated. We first give the expression for Pr(T N,k ≤ t), the probability that at least k links are generated before t. This is equivalent to the probability that, at time t, the number of elementary links that have not yet been generated is N − k or less 149 : where Pr(T 0 ≤ t) = 1 − (1 − p g ) t since T 0 is geometrically distributed with success probability p g . The probability that precisely k of the N segments are generated at time t is from which the mean waiting time is calculated as The mean waiting time T N,k from Eq. (62) can be computed by solving a recurrence formula where T N,k is expressed as function of T N,k−1 150,151 . For k = N, i.e. the waiting time that all elementary entanglement are prepared 150 , the solution reads For p g 1, this expression can be approximated by 142,144,149 with where H(N) is the N-th harmonic number and γ ≈ 0.57721 is the Euler-Mascheroni constant. In separate work, Praxmeyer included finite memory time with cut-off into the calculation and obtained 151 where τ is the cut-off threshold and q g = 1 − p g . Similar derivations as the ones above can be used for the waiting time until the first, instead of the last, elementary link has been generated. Those derivations are relevant for the analysis of multiplexed repeater protocols and the mean waiting time in those cases has been analyzed in 125,149 .
To our knowledge, in contrast to the waiting time, there is no exact fidelity calculation with exponential memory decoherence for deterministic swap. A lower bound on the fidelity can be obtained by assuming the worst case, i.e. the swap is performed only after all elementary links are generated 149,150 .
These expressions presented here for the deterministicswap case also apply to repeater chains where the numbers of segments is not a power of 2, as well as to more general network topology 152 . The reason for this is that if the swaps are deterministic, the waiting time equals the time until all elementary links in the network have been prepared. Thus, the nested structure does not exist and the only relevant parameters are the number of elementary links and the elementarylink success probability p g .
The waiting time in the deterministic-swap case can be used as an approximation to the case where p s is slightly lower than 1 and bounds from below the waiting time for general p s . The quality of the approximation is shown in Fig. 4.

Methods for analyzing distillation-based repeater schemes with memory-decoherence
In contrast to entanglement swapping, distillation has a fidelity-dependent success probability. In the absence of decoherence, the fidelity of a pair does not decrease while it is waiting for other components to succeed. Hence, the success probability p d is a constant for each level and distillation can be studied in the same way as entanglement swapping. However, in the presence of decoherence, fidelity and success probability become correlated, which complicates the analysis.
We finish by mentioning two tools for bounding the fidelity and generation rate of distillation-based repeater schemes in the presence of memory decoherence. First, upper bounds on the achievable fidelity can be derived using fixed-point analysis 122,127 . In this approach, one makes use of the fact that entanglement distillation does not improve the fidelity when the quality of the input links is sufficiently high. Such a fidelity is thus a fixed point of the entanglement distillation procedure and it depends on the quality of the local operations 127 and memories 141 . If the fixed-point is an attractor and the input links have fidelity lower than the fixed-point, repeated application of entanglement distillation cannot boost fidelity beyond the fixed-point and it thus forms an upper bound. Next, lower bounds on the fidelity can be trivially obtained for protocols that impose a fidelity cut-off, i.e. protocols that discard the entanglement if the fidelity is lower than a certain threshold 153 . Because the distillation success probability is a monotonic function of the fidelity of the input states, a lower bound on the fidelity by a cut-off also directly yields a lower bound on the success probability.

C. Numerical tools for evaluating analytical expressions
Above, in Section III B, we reviewed analytical tools for computing the probability distribution of the waiting time for generating remote entanglement and of the entanglement quality. For models that do not include memory decoherence, distillation, or cutoffs, these tools are sufficient. For more complex models and for the analysis of many-node networks, the tools presented above are often still applicable but analytically evaluating the resulting expressions to compact, closedform expressions is unfeasible. An example of such a case is a nested repeater chain with cut-offs and non-deterministic swapping. In this case, no concise analytical expression for the waiting time is known. A priory, it is possible to write down a recursive analytical expression for the waiting time using a similar reasoning to Section III B 2, where the singlerepeater case was treated. However, the recursion relation has so far not been solved for general repeater chains. Fortunately, numerical tools enable the evaluation of such expressions. In this section, we treat three classes of such tools: Markov Chain algorithms, probability-tracking algorithms and Monte Carlo methods for abstract models. These numerical tools differ from the simulation tools in Section IV in that they are all directly applied to expressions that can in principle be evaluated analytically.

Markov Chain methods
In many repeater protocols, the change of the entanglement in the network in the next time step only depends on the existing entanglement. Shchukin et al. used this idea to model the network as a discrete Markov chain 144 , which can be visually depicted as a directed graph, an example of which is shown in Fig. 5. A vertex in this graph is a state of the network, i.e. the collection of entanglement that exists at a given point in time. The network transitions from one state to the other with a fixed probability, which is visualized by directed edges of the labeled graph. At each time step, a network randomly transitions from its current state to the next state according to the transition probabilities over outgoing edges. For example, in the single-repeater protocol depicted in Fig. 5, the transition from the 'an entangled pair exists on each of the two segments' state (11) to 'entanglement exists between end nodes' (11) occurs with the entanglement swapping success probability p s (entanglement swapping succeeded), whereas the transition to 'no entanglement' (state 00) has probability 1 − p s (entanglement swapping failed and the two involved links are lost).
An equivalent representation of a Markov chain is the transition probability matrix (TPM), where entry ( j, k) represents the transition probability from state j to state k. Since a single transition corresponds to a single time step, the waiting time distribution equals the distribution of the number of edges traversed before reaching a predefined target state, such as 'entanglement between the end nodes of the repeater chain' (11 in  (1 − p g ) 2 The original proposal by Shchukin et al. included an analysis of the waiting time. Later, the idea was refined to include memory decoherence by Vinay et al. 153 They computed the fidelity distribution by assigning a noise parameter to certain transition edges and calculated how many times the edges are traversed given that the entanglement distribution is completed at time t. With this noise model, the Markov chain method was used to study the BDCZ protocol 127 , which includes entanglement distillation. Due to the assumption of the Markov process, i.e. the system has no memory of the past, this method cannot handle fidelity-dependent success probability without assigning each possible fidelity a state representation. As an alternative, Vinay et al. provided a lower bound to the final fidelity using fidelity cutoffs (see Section III B 5).
Apart from repeater chain protocols, Markov chains have also been used to study more general network protocols, such as a quantum switch by Vardoyan et al. 154,155 , who used the continuous-time Markov chains as an approximation to discrete-time Markov chains. In this model, the transition probability is replaced by the transition rate. Compared to their discrete counterparts, continuous-time Markov chains simplify the analysis in various aspects. For instance, Vardoyan et al. included a model of decoherence where the states are lost at a fixed rate by adding an additional transition edge indicating the loss of one entangled pair. Moreover, Vardoyan et al. show that the quality of the approximation is high in many scenarios 154 .
The Markov chain method is rather general and flexible: in principle, the waiting time of any entanglement distribution protocol with predefined transition probabilities can be calculated, regardless of the topology or entanglement swapping policy (such as swapping as soon as two links are available, regardless of the segments on which this entanglement has been produced, or swapping only between predefined segments). However, this method is computationally expensive. The size of the TPM is the same as the number of possible Markov states and, in general, grows exponentially with the number of nodes.
This rapid growth of the size of TPM can be partially mitigated. For instance, by grouping equivalent Markov states and treating them as one state, the complexity can be drastically Finally, in a recent development, Khatri 156 introduced a method for describing network protocols based on quantum partially observable Markov decision processes 157 . A quantum partially observable Markov decision process is a reinforcement-learning based framework for protocol optimization. In this framework, the protocol obtains feedback from its actions in the form of classical information about the quantum state that the network holds, which it uses to optimize the next action it will perform. As an application of the method, Khatri found analytical solutions for optimizing a cut-off for elementary link generation under different constraints. It is an interesting open question whether this approach can be extended for efficiently characterizing and optimizing protocols over large repeater chains and networks.

Probability-tracking algorithm
The next numerical tool tracks the full waiting time probability distribution and the average fidelity of the delivered quantum state. We explain this method via a concrete example, a symmetric nested repeater protocol with 2 n segments and no entanglement distillation (depicted in Fig. 2(b) for n = 2). In Section III B 2, we treated the case for n = 1, which resulted in an expression for the waiting time random variable consisting of the maximum of two copies of the waiting time of the bottom level (Eq. (51)) and a compound sum (Eq. (54)). The first element, the maximum, corresponds to the fact that an entanglement swap acts on two links that need both be generated, so one needs to wait until the latest of the two links has been prepared. The second element is the sum of the waiting time until the first successful swap attempt. Since the number of attempts is probabilistic, the result is a compound sum. The analysis for the n = 1 case can be generalized to an arbitrary number of nesting levels n and yields a recursive expression of the waiting time T n which alternates between compound sums and maximums of two copies of the waiting time T n−1 on the previous level. Unfortunately, as discussed in Section III B, to our knowledge, this recursive expression of T n has not been analytically evaluated for n > 1. Hence, various approximation methods were introduced in that section. The exact evaluation can, however, be achieved with numerical tools. By truncating the waiting time at a prespecified time t trunc , the waiting time distribution becomes finite. The evaluation with numerical tools leads to an algorithm that tracks both the truncated probability distribution of T n and the associated fidelity 158,159 .
On the 2 n -segment nested repeater protocol, the algorithm computes the waiting time distribution as follows. If n = 0, i.e. if there is no repeater and the two end nodes obtain entanglement by direct generation, the waiting time T 0 follows a geometric distribution (Eq. (49)). If n = 1, i.e. there is a single repeater and two segments, then the algorithm evaluates Eq. (54), which describes how the probability distribution of the waiting time T 1 can be obtained from the distribution of T 0 . Although the two elements in Eq. (54), the maximum and the geometric compound sum, can in principle be evaluated sequentially 158 , a computationally faster approach is to separate the probability distributions of failed and successful swap attempts 159 . For n > 1, the algorithm is applied iteratively over the nesting levels until level n has been reached. The algorithm can be extended in polynomial time to also track the average fidelity of the produced quantum state. This fidelity is a function of the delivery time and it can include the effect of memory decoherence.
Although the example protocol above only consists of entanglement swapping, the algorithm is applicable to any protocol which is composed of entanglement distillation and cutoffs, in addition to entanglement swaps 159 . The algorithm presupposes that the protocol is composed of these three operations in a predefined order, e.g. which entangled pairs are swapped must be known in advance. The algorithm scales polynomially in the number of nodes and in the truncation time t trunc and has been used to track over 1000 nodes for some parameter regimes 158 .
A related approach to the probability-tracking algorithm explained above is taken by Kuzmin et al. 145,146 . This method assumes that the waiting time of an elementary link is exponentially distributed (Eq. (50)), after which the mean waiting time for the next level is computed by evaluating a continuous integral, as well as the quantum state in the presence of mem-ory decay. These are then used as input to the next nesting level, by assuming that at that level, the waiting time follows an exponential distribution also. With this approximation, the calculation of the maximum in Eq. (51) is simplified.

Sampling the analytical expressions with Monte Carlo method
So far in this section, we have discussed two methods that compute the statistical distribution of the waiting time and produced quantum state. For a given model, both of them evaluate the analytical expression exactly up to machine precision. An alternative to this semi-analytical computation is to sample the expressions on random variables for the waiting time and the delivered state using a Monte Carlo approach 158 . Instead of tracking the whole distribution, this approach samples a pair of waiting time and the produced state between the end nodes. By sampling many times, the probability distribution of the waiting time and the quantum state can be reconstructed.
Again, let us take the 2 n -segment nested repeater protocol as an example to explain the algorithm. The individual sample pairs are produced by iterating over the different components of the repeater protocol, following its nested structure. At each component, a pair is sampled recursively, following the expressions on random variables, which thus become expressions on individual events. For instance, Eq. (51) requires sampling two instances of T 0 for entanglement generation and then taking the maximum of both to produce a sample of M 0 . Also, memory decoherence can be calculated from the time difference of two events. Similarly, the method can handle cut-offs. Furthermore, distillation can also be included in the protocol, since the input states to a distillation attempt, which determine its success probability, are also sampled. For nested protocols, the Monte Carlo algorithm can be defined as a recursive function, following the nested structure of the protocols.

D. Second and third generation repeater protocols
So far, we have only treated first-generation quantum repeater protocols, i.e. protocols for which the building blocks -fresh entanglement generation, entanglement swapping, and entanglement distillation -are probabilistic operations. The quantum repeater proposals that do not fall into this category make explicit use of quantum error correction codes and are referred to as second-generation (probabilistic entanglement generation, deterministic entanglement swapping, and oneway quantum error correction) and third-generation repeaters (loss-tolerant entanglement generation) 114,160 .
In first-generation repeaters, entanglement generation and swapping are probabilistic, and once it has succeeded, the entanglement is kept in quantum memories and needs to wait until other components performed in parallel have succeeded. Consequently, the waiting time probability distribution and state quality are a complex function of the success probabilities of components in the repeater chain.
In contrast, for second-generation repeater protocols, such as [161][162][163][164] , entanglement swapping and (one-way) quantum error correction are no longer probabilistic (although entanglement generation is still probabilistic, it may be parallelized to achieve near-unit generation success probability). As a result, the time at which the entire repeater chain finishes with a single attempt at generating end-to-end entanglement is simply a sum of the (constant) times that the individual components take. Not only there is no waiting for other components, but there are also no feedback loops here, i.e. components that need to restart in case others have failed. The unit time step at which such repeater chains operate is an attempt at end-to-end entanglement generation (i.e. the sum of the individual component times); the probability that such an attempt succeeds is the product of all individual steps succeeding. For this reason, the distribution of the waiting time is a geometric distribution.
A similar reasoning applies to third generation repeaters 160,[165][166][167][168][169][170][171][172][173][174][175] , where entanglement between adjacent nodes is established almost deterministically, rather than probabilistically, by encoding part of locally-generated entanglement into a large state of photons, followed by transmission of the encoded state. Commonly, the analysis of the propagation of operational errors (for 2nd and 3rd generation) and the propagation of physical loss errors into logical errors (for 3rd generation) is based on work on quantum error correction codes (combined with explicit counting of the combinations of losses of photons which yield errors beyond recovery) and optical quantum computation. We consider such tools out of scope for this review.

IV. SIMULATION TOOLS
In this section, we sketch the methodology for evaluating the network performance based on simulation and then discuss the particularities of quantum networks together with a selection of existing tools for quantum networks. For an indepth discussion about network simulation, we refer to the books [176][177][178] .
Let us first compare simulation tools against analytical tools. Analytical tools are well suited for predicting the performance of scenarios which are simple in terms of streamlined models, but also in terms of the communication protocol, network topology, and network usage. These scenarios typically are simplified versions of the scenarios of interest and the output of analytical tools can be understood as a rough statement about feasibility (Section III) or of the optimal achievable performance (Section II). Simulation finds its role in use cases that move beyond simple scenarios.

A. Simulation methodology
In the following, we sketch a methodology for network performance evaluation via simulation. The methodology consists of four steps: problem formulation, theoretical modeling, implementation, and simulation. It is a streamlined version of the best practice methodologies discussed in 1,176 .
Problem formulation. The first step is to define the goal of the performance evaluation. Some examples are the validation of protocol designs, the comparison between several alternative solutions, the study of network stability, the sensitivity analysis of performance with respect to hardware parameters, parameter optimization, etc. After the problem has been formulated, the desired performance metrics should be identified; different metrics impose different constraints to the modeling step. Some examples are: average rate, fidelity, latency, throughput, etc.
Theoretical modeling. The second step is modeling. Best practices 1,176 dictate a separation between the theoretical characterization of a network element and its implementation. The theoretical characterization begins with an identification of boundary conditions and assumptions. An important design consideration is the level of detail of the models and protocols. This is a difficult choice that depends on several factors in addition to the problem goal and the performance metric: the availability or lack of data can condition the modeling of some elements, and the computation time might limit the level of detail as well as more general constraints.
Let us sketch the qualitative trade-offs between simple and detailed modeling. Simple models have several advantages: the modeling effort is reduced, simulations are comparatively faster and they can be easily modified. Conversely, the output should be taken as an estimation of the real behavior of the network. In some cases, the output of the simulations might be similar to that of analytical methods. As the models increase in detail, the simulations become slower, it becomes more difficult to modify the network and the output becomes representative of the output expected from the real devices. As the level of detail increases, the cost also increases.
The output of the theoretical modeling process is a full specification of the behavior of all elements of interest together with their interactions. The relevant question at this stage is whether the level of detail is both necessary and sufficient. That is, whether either the model does not allow to evaluate the desired metric or, on the contrary, a simpler model would be sufficient. This can be quantified by estimating whether additional or smaller amounts of detail will impact the simulation outcome. Finally, the model is documented.
Implementation. The documented model is then translated into software. The first decision is whether to develop an ad-hoc solution or to opt for an existing framework. The advantage of an ad-hoc solution is that it can be tailored to the project at hand, which in turn can benefit both the accuracy of the implementation and its computational efficiency. On the other hand, tailored solutions have several drawbacks: they might be difficult to expand at a later stage, they can be costly to develop and it is difficult to compare the outcome with other simulation studies. Most projects build on an existing framework, and the choice can depend on many factors including the existence of models, the validation means, the capability of quantifying the desired metric, the cost, the flexibility for further investigation, and the accuracy.
Simulation accuracy has a particular implication for the  choice of the simulation engine. If accurate physical modeling is a requirement, accurate characterization of timing between actions is necessary. Because of this reason, it is already the case that most classical network simulators are based on the discrete-event paradigm 176 , a paradigm for simulation that guarantees reliable timing behavior (see Figure 6). The requirement of accurate tracking of time is particularly strong in the case of quantum networks, for instance, if an accurate characterization of the average fidelity is a goal; the reason is that the amount of noise a quantum memory undergoes is highly dependent on time. In consequence, inaccurate time handling can yield strong variations in the performance evaluation. A side benefit of discrete event simulators is that simulations are inherently repeatable for a fixed randomness seed. This is more difficult for simulators where the ordering of events is not guaranteed, such as in multi-threaded or distributed implementations.
Once the models have been translated to software, they should first be verified against the theoretical models and, if possible, validated against experimental data.
Simulation. The final step is the performance of the simulations and the analysis of the results.

B. The challenge of simulating quantum networks
The main difference between simulating quantum and classical networks is that in full generality the state of the network can be an entangled state. This implies that the state can not be described as an aggregate of the local description of the states at each network node. Moreover, the representation of the state grows exponentially with the number of qubits in the network. For this reason, classical simulators of quantum computation can not simulate computations involving more than tens of qubits 179,180 .
Fortunately, many scenarios of interest in quantum networks involve only Clifford gates and Pauli noise, yielding classical simulation tractable 181 . Other scenarios involve many small entangled states which dynamically interact, shrink via measurement or are fused but never grow larger than a few qubits. For these scenarios, classical simulation is also possible.

C. Existing simulation platforms
We now introduce several existing platforms for simulating quantum networks. We limit our discussion to platforms that can simulate arbitrary quantum network topologies and protocols. There are a number of specialized platforms that fall beyond the scope of this review [182][183][184][185][186][187][188][189][190][191] .
The first two platforms are SimulaQron 192 and the QUantum NETwork SIMulator (QuNetSim) 193 . Both of them aim at facilitating quantum network application development. QuNetSim focuses on ease of implementation. It posits a concrete architecture of quantum networks inspired by the OSI model 194 and a convenient feature of the architecture is the joint arrival of the quantum payload together with a corresponding classical header. This simplifies the process of synchronization and the logic at the nodes, because a node's action after message reception can now be based on the classical header, irrespective of whether the payload is classical or quantum. QuNetSim has been used to explore routing schemes 195 . SimulaQron can be run distributed over a classical network, i.e. on physically-distinct machines, to simulate a network of quantum computers. The back-end can be accessed via a quantum hardware interface, thus acting as a substitute for the quantum device. SimulaQron has been used to explore entanglement verification schemes 196 and several multipartite protocols 197,198 ; moreover, within the context of the quantum protocol zoo 199 , SimulaQron has been used to implement other protocols such as a quantum cheque scheme 200 , a protocol for leader election 201 , universal blind quantum computation 202 and coin flipping 203 .
The next platform is the Simulator for Quantum Networks and CHannels (SQUANCH) 204 . The goal of SQUANCH is to simulate quantum information protocols with realistic noise models over large networks. The simulator builds on the notion of the agent, a party in the quantum network connected to other agents via classical and quantum channels. An agent can manipulate in runtime a classical and a quantum memory. The quantum memory can be a subset of a larger distributed quantum state, and the underlying registers are dynamically resized as they get measured or fused. The simulation mirrors the network configuration by assigning a different process to each agent, and performance is optimized with a NumPy backend. It has been used to evaluate quantum channel discrimination schemes 205 and to investigate a photonic lattice architecture for a universal quantum programmable gate array 206 .
While noise can be added to these first three simulators, and noise itself can sometimes be accurately modeled, the noise parameter depends on the time elapsed which is challenging to capture. Next, we introduce three discrete event simulators.
The QUantum Internet Simulator Package (QuISP) 207,208 aims to capture complex quantum network behavior. In particular, one of the target use cases is understanding emergent behavior with complex networks composed of large quantum networks interconnected. QuISP has been designed as a module for OMNeT++, an open-source simulator of classical networks with a long history 209 . This implies that most of OM-NeT++ functionality is potentially available. This includes visualization, logging, and data analysis tools but also the models and protocols of classical networks. The main mode of operation of QuISP is based on tracking Pauli errors. While this is not accurate enough for link-layer protocols, it enables very efficient simulation. QuISP has already been used to study the capability of first-generation repeater links including realistic parameters compatible with Pauli tracking and to evaluate the RuleSet-based communication protocol 13 .
The two last discrete event simulators are the Simulator of QUantum Network Communication (SeQUeNCE) 210 and the NETwork Simulator for QUantum Information with Discrete events (NetSquid) 211 . Both projects aim at accurate physical simulation. SeQUenCE is designed according to a concrete network stack architecture, consisting of reprogrammable modules for given functionalities. SeQUeNCE has been used for a full simulation of the BB84 quantum key distribution protocol and quantum teleportation 212,213 , as well as for analysis of a nine-node quantum network with quantum memories modeled after single erbium ions in solids 210 . Regarding NetSquid, some particular features are a dynamical handling of the quantum registers similar to SQUANCH and a fully parallelizable architecture. NetSquid has shown to enable physically-accurate modeling through simulation of a link layer protocol with node models based on NV centers in diamond 8 . Its modular design was showcased by comparing two types of atomic-ensemble based quantum memories through only replacing the hardware model component, while its quantum engine is powerful enough to simulate up to one thousand nodes for some models 211 . NetSquid has also been used to evaluate a quantum router architecture 214 and to go beyond the analytical regime in the analysis of a quantum switch 154,211 .

V. OUTLOOK
Quantum networks will enable the implementation of communication tasks with qualitative advantages with respect to the communication networks we know today. To compare different solutions, optimize over parameter space and inform experiments, it is necessary to evaluate the performance of concrete quantum network scenarios. Here, we discussed the existing tools for evaluating the performance of quantum networks from three different angles: information-theoretic benchmarks, analytical tools, and simulation.
First, we presented information-theoretic tools, which allow comparing the performance that a protocol with realistic parameters achieves against the optimal performance of a protocol with similar resources. In particular, we discussed a number of recent tools that allow us to bound the capacity of the network for entanglement and key distribution in multi-user scenarios with a complexity that grows polynomially with the number of nodes in the network.
Then we discussed analytical tools, mostly for estimating the average waiting time and fidelity of entanglement distribution over repeater chains. We described tools that can analyze moderately complicated protocols composed of entanglement generation, distillation, swapping, and cut-offs, and discussed the trade-offs between different approximations used in the literature.
Simulation tools take the limelight when one moves away from simple scenarios. They can accurately model complicated scenarios and be used for problems such as validate protocols, evaluate alternative network solutions, study the stability of a network, etc. Here, we reviewed the methodology for network simulation in the context of quantum networks and introduced six simulation frameworks with different goals and target audiences.
Our goal in this review has been to cover the tools for evaluating the performance of quantum networks. This is a rapidly developing topic. For instance, most of the benchmarks for quantum networks have been developed in the past five years and all the six simulation frameworks discussed have been released or presented to the scientific community in the past two years. We hope that this review raises awareness of the wealth of existing tools and stimulates the development of new ones.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding authors upon request.