Quantum networks will enable the implementation of communication tasks with qualitative advantages with respect to the communication networks known 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, the authors review the state-of-the-art of tools for evaluating the performance of quantum networks. The authors present them from three different angles: information-theoretic benchmarks, analytical tools, and simulation.

## 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 a 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 cannot be copied, reducing the applicability of classical approaches for long-distance communication. 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 Sec. 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 communication 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 nonvanishing 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, then this certifies that the device qualitatively improves the communication capabilities of the two parties.

Second, in Sec. III, we review analytical tools for evaluating the performance of concrete protocols. We focus our analysis on 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 cutoffs.

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 Sec. 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 Sec. 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 and optimization of quantum networks^{7–20} beyond the scope of this review. We also omit the closely related topics of quantum networks for creating correlations^{21–25} and of complex quantum networks.^{26,27} Finally, there is a large literature regarding communication 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.^{28–30} 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.^{31,32} Quantum channels can also be abstracted by a mathematical object (see below Sec. 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.^{33} 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 multiuser 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.^{30,34} 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 state $\rho \u0302$ as an input and returns a quantum state $\sigma \u0302k$ with probability *p _{k}* as an output, where $\sigma \u0302k:=M\u0302kX\rho \u0302(M\u0302kX)\u2020/pk$ with Kraus operators ${M\u0302kX}k$ satisfying $\u2211k(M\u0302kX)\u2020M\u0302kX=1\u0302X$ and $pk=Tr[(M\u0302kX)\u2020M\u0302kX\rho \u0302]$. 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 fiber, a superconducting microwave transmission line, or an optical free space link. In theory, any quantum channel $nX\u2192Y$ from a node

*X*to a node

*Y*is described by a completely positive trace-preserving (CPTP) map. In particular, quantum channel $nX\u2192Y$ for an input state $\rho \u0302$ is described by $nX\u2192Y(\rho \u0302):=TrE\u2032[U\u0302XE\u2192YE\u2032(\rho \u0302\u2297|0\u27e9\u27e80|E)(U\u0302XE\u2192YE\u2032)\u2020]$ with a unitary operator $U\u0302XE\u2192YE\u2032$ on Hilbert space $hX\u2297hE=hY\u2297hE\u2032$ and a state $|0\u27e9E$ 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. 1(a)]. A quantum information processing node in a quantum network is regarded as a vertex *v* in a graph, while a quantum channel $nv\u2192v\u2032$ to send a quantum system from a node *v* to another node $v\u2032$ is associated with a directed edge $e=v\u2192v\u2032$ [or, simply, $e=vv\u2032$, where the former *v* and latter $v\u2032$ 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(\u2208V)$ wants to send an unknown quantum state with dimension *d* to another client $B(\u2208V)$ via quantum teleportation,^{36} they need to share a bipartite maximally entangled state $|\Phi d\u27e9AB:=\u2211i=1d|ii\u27e9AB/d$. Even if the client *A* wants to send bits to the client *B* secretly, since the maximally entangled state provides them $log2d$ bits of secret key for the one-time pad, sharing the maximally entangled state is sufficient. Indeed, there is a quantum key distribution (QKD) scheme based on sharing the maximally entangled state. However, in general, a maximally entangled state is not necessary for QKD. In fact, there exists a larger class of so-called private states^{37,38} $\gamma \u0302dAB$, from which clients *AB* can obtain $log2d$ secret bits. Such states are of the form $\gamma \u0302dAB:=U\u0302twistAB(|\Phi d\u27e9\u27e8\Phi d|A\u2032B\u2032\u2297\sigma \u0302A\u2033B\u2033)(U\u0302twistAB)\u2020$, with an arbitrary shared state $\sigma \u0302A\u2033B\u2033$ and an arbitrary controlled unitary operation in the form $U\u0302twistAB:=\u2211i,j=1d|ij\u27e9\u27e8ij|A\u2032B\u2032\u2297U\u0302ijA\u2033B\u2033$. More generally, one could think of a group of clients $C:=C1C2\cdots Cn(\u2282V)$ ($C1,C2,\u2026,Cn\u2208V$) who want to share a multipartite entangled state, such as a Greenberger–Horne–Zeilinger (GHZ) state^{39} $|GHZd\u27e9C:=\u2211i=1d|ii\cdots i\u27e9C1C2\cdots Cn/d$, which can be used for secret sharing,^{40} a clock synchronization protocol^{4} as well as quantum conference key agreement.^{41} In the case of quantum conference key agreement, there is, again, a larger class of so-called multipartite privates $\gamma \u0302dC$, from which the group of clients can obtain $log2\u2009d$ bits of conference key.^{41} Such states are of the form $\gamma \u0302dC:=U\u0302twistC1C2\cdots Cn(|GHZd\u27e9\u27e8GHZd|C1\u2032C2\u2032\cdots Cn\u2032\u2297\sigma \u0302C1\u2033C2\u2033\cdots Cn\u2033)(U\u0302twistC1C2\cdots Cn)\u2020$, where $U\u0302twistC1C2\cdots Cn$ is an arbitrary controlled unitary operation in the form $U\u0302twistC1C2\cdots Cn:=\u2211i1,i2,\u2026,in=1d|i1i2\cdots in\u27e9\u27e8i1i2\cdots in|C1\u2032C2\u2032\cdots Cn\u2032\u2297U\u0302i1i2\cdots inC1\u2033C2\u2033\cdots Cn\u2033$ and $\sigma \u0302C1\u2033C2\u2033\cdots Cn\u2033$ is an arbitrary shared state.

Depending on the communication task that clients want to perform, clients *C* may ask a provider to supply a resource entangled state $\tau \u0302dC$ for the task. For example, $\tau \u0302dC$ may be a bipartite maximally entangled state $|\Phi d\u27e9C$, a private state $\gamma \u0302dAB$, a GHZ state $|GHZd\u27e9C$, or a multipartite private state $\gamma \u0302dC$. The amount of the entangled resource present in $\tau \u0302dC$ is quantified by the parameter *d*. Namely, $log2d$ 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 Sec. 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\u2282V$. This is because entanglement cannot be generated only with LOCC in general (or, precisely, entanglement is defined^{42} as a state which cannot be generated by LOCC, that is, which cannot be described by the form of fully separable states $\u2211iqi\u2297v\u2208V\rho \u0302iv$, where ${qi}i$ is a probability distribution and $\rho \u0302iv$ is a quantum state of node $v\u2208V$), and thus, the provider needs to use a quantum network. Without loss of generality, a protocol of the provider which provides a final state $\rho \u0302klV$ close to a target state $\tau \u0302dklC$ with probability $pkl$ is described as follows [see also Fig. 1(b)]. A provider determines an available quantum network, associated with a graph $G=(V,E)$ such that $C\u2282V$. The protocol starts by preparing the whole system in a fully separable state $\rho \u03021V$ and then by using a quantum channel $ne1$ with $e1\u2208E$. This is followed by arbitrary LOCC among all the nodes, which gives an outcome *k*_{1} and a quantum state $\rho \u0302k1V$ with probability $pk1$ (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 $nek1$ with $ek1\u2208E$, followed by LOCC among all the nodes. This LOCC gives an outcome *k*_{2} and a quantum state $\rho \u0302k2k1V$ with probability $pk2|k1$. Similarly, in the *i*th round, according to the previous outcomes $ki\u22121:=ki\u22121\u2026k2k1$ (with $k0:=1$), the protocol may use a quantum channel $neki\u22121$ with $eki\u22121\u2208E$, followed by LOCC providing a quantum state $\rho \u0302kiV$ with a new outcome *k _{i}* with probability $pki|ki\u22121$. Eventually, after a finite number of rounds, say after an

*l*-th round, the protocol presents clients

*C*with $\rho \u0302klC=TrV\u2216C(\rho \u0302klV)$ close to a target state $\tau \u0302dklC$ in the sense of $||\rho \u0302klC\u2212\tau \u0302dklC||1\u2264\epsilon $ for $\epsilon >0$, with probability $pkl:=pkl|kl\u22121\u2026pk3|k2pk2|k1pk1$, where $||X\u0302||1$ is the trace norm defined by $||X\u0302||1:=TrX\u0302\u2020X\u0302$.

Important quantities related to the evaluation of a protocol are (1) how many times quantum channels $ne$ are used in the protocol, (2) how valuable the state $\tau \u0302dC$ 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\xafe:=\u2211i=0l\u22121\u27e8\delta e,eki\u27e9ki=\u2211i=0l\u22121\u2211kipki\delta e,eki$ with the Kronecker delta $\delta i,j$ and $l\xafE:=\u2211e\u2208El\xafe$ are important quantities because $l\xafe$ represents the average number of times quantum channel $ne$ is used in the protocol and $l\xafE$ represents the average number of total channel uses. Thus, ${l\xafe}e\u2208E$ can be a set to characterize the protocol. If we introduce an average usage frequency $f\xafe:=l\xafe/l\xafE(\u22650)$ satisfying $\u2211e\u2208Ef\xafe=1$, a set ${l\xafE,{f\xafe}e\u2208E}$ can be used to characterize the protocol. Or, by introducing an average usage rate $r\xafe:=l\xafe/l(\u22650)$, a set ${l,{r\xafe}e\u2208E}$ can also be used to characterize the protocol, where *l* is a parameter like time (because it has no explicit relation with $l\xafe$ in contrast to $l\xafE$). As for (2), in general, it is not so obvious to define a quantity for evaluating the value of the target state $\tau \u0302dklC$ because the value depends on the final application in general and may be determined by a complex function of the target state $\tau \u0302dC$ or its parameter *d*. However, if the target state $\tau \u0302dC$ is GHZ state $|GHZd\u27e9C$ or private state $\gamma dC$, the value of the target state can be quantified by $log2d$ 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 $\tau \u0302klC$ the final state $\rho \u0302klC$ is. In general, this trace distance $||\rho \u0302\u2212\sigma \u0302||1$ is related^{43} with another well known measure, the fidelity $F(\rho \u0302,\sigma \u0302):=||\rho \u0302\sigma \u0302||1$, as $2(1\u2212F)\u2264||\rho \u0302\u2212\sigma \u0302||1\u226421\u2212F$. 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,^{44} in contrast to the fidelity.

#### 1. 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^{45,46} and Azuma *et al.*^{35} In particular, in this case, the target state $\tau \u0302dC$ is regarded as a single bipartite state $|\Phi d\u27e9AB$ or $\gamma \u0302dAB$, where $A,B\u2208V$. To obtain such upper bounds, Pirandola generalizes the Pirandola–Laurenza–Ottaviani–Banchi (PLOB) bound^{33} on the private capacity for point-to-point quantum key distribution with a pure-loss bosonic channel, while Azuma *et al.* use the Takeoka–Guha–Wilde (TGW) bound^{47} for that. The shared methodology is summarized by Rigovacca *et al.*,^{48} which is described in what follows.

An entanglement measure^{42} $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: (i) If a state $\rho \u0302XY$ is close to a target state $\tau \u0302dXY$ which is $|\Phi d\u27e9\u27e8\Phi d|XY$ (or $\gamma \u0302dXY$), that is, $||\rho \u0302XY\u2212\tau \u0302dXY||1\u2264\epsilon $, then there exist two real continuous functions $fE(\epsilon )$ and $gE(\epsilon )$ with $lim\epsilon \u21920fE(\epsilon )=0$ and $lim\epsilon \u21920gE(\epsilon )=1$, such that $E(\rho \u0302XY)\u2265gE(\epsilon )\u2009log2d\u2212fE(\epsilon )$; (ii) we have $E(nX\u2032\u2192Y\u2032(\rho \u0302XX\u2032Y))\u2264E(nX\u2032\u2192Y\u2032)+E(\rho \u0302XX\u2032Y)$ for any state $\rho \u0302XX\u2032Y$, where $E(nX\u2192Y):=max\sigma \u0302XX\u2032E(nX\u2192Y(\sigma \u0302XX\u2032))$ is called the entanglement of channel.^{33,47,49,50} 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(\u2208V)$ and $B(\u2208V)$ with maximally entangled state $|\Phi dkl\u27e9AB$ (or private state $\gamma \u0302dklAB$) with probability $pkl$ and error $\epsilon >0$, by using a quantum network associated with graph $G=(V,E)$, it holds

for any $VA\u2282V$ which is a set of nodes with $A\u2208VA$ and $B\u2209VA$, where $\u27e8xk\u27e9k:=\u2211kpkxk$ and $\u2202(VA)(\u2282E)$ is the set of edges which connect a node in *V _{A}* and a node in $V\u2216VA$ (Fig. 1). Since Eq. (1) holds for any choice of

*V*, the minimization of the right-hand side of Eq. (1) over the choice of

_{A}*V*gives the tightest bound on the left-hand side. In general, there might be cases where we have no computable expressions for point-to-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 use $l\xafe$ of each channel $ne$. In contrast, the upper bound (1) is always estimable, once we are given ${E(ne)}e\u2208E$ or their upper bounds, because it depends only on a single use of channel $ne$, rather than multiple uses of it.

_{A}Examples of entanglement measure $E$ satisfying property (i)^{48,49,51} for either $|\Phi d\u27e9AB$ or $\gamma \u0302dAB$ and property (ii)^{47,49,52} are the squashed entanglement^{53,54} $Esq$ and the max-relative entropy^{55} $Emax$ of entanglement. The relative entropy^{56,57} $ER$ of entanglement satisfies^{33,50,58} property (i) for either $|\Phi d\u27e9AB$ or $\gamma \u0302dAB$ in general, but, at present, it is shown^{33,45,50} to satisfy property (ii) only for some quantum channels, called teleportation simulable/stretchable channels or Choi simulable channels.^{33,59–65} A channel $nA\u2192B$ is called teleportation simulable if there is a deterministic LOCC operation $lA\u2032A\u2033:B$ between systems $A\u2032A\u2033$ and *B* such that

for any state $\rho \u0302A$, where $d=dimhA$. Examples of these channels are ones which “commute” with the unitary correction in the quantum teleportation. In particular, the quantum teleportation^{36} from system $A\u2033$ in an unknown state $\rho \u0302A\u2033$ to system *A* is implemented by generalized Bell measurement ${M\u0302iA\u2032A\u2033}i=1,2,\u2026,d2$ on systems $A\u2032A\u2033$ in state $|\Phi d\u27e9\u27e8\Phi d|AA\u2032\u2297\rho \u0302A\u2033$, followed by a unitary operation $UiA$ depending on the random outcome *i* of the Bell measurement. This means $\rho \u0302A=d2(U\u0302iA\u2297M\u0302iA\u2032A\u2033)(|\Phi d\u27e9\u27e8\Phi d|AA\u2032\u2297\rho \u0302A\u2033)(U\u0302iA\u2297M\u0302iA\u2032A\u2033)\u2020$. Therefore, if every correction unitary $U\u0302iA$ “commutes” with $nA\u2192B$, more precisely, if there is a unitary operation $V\u0302iB$ with $nA\u2192B(U\u0302iA\sigma \u0302A(U\u0302iA)\u2020)=V\u0302iBnA\u2192B(\sigma \u0302A)(V\u0302iB)\u2020$ for any state $\sigma \u0302A$ and any outcome *i*, then the channel $nA\u2192B$ is teleportation simulable. For example, depolarizing channels, phase-flip 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\u2192B}$), the definitions of quantum/private capacities in a quantum network themselves may have variety.^{66,67}

A rather general definition is

for given ${qe}e\u2208E$ with $qe\u22650$, where $n\u22650,\u2009\Lambda (n,{qe}e\u2208E,\epsilon )$ is a set of protocols which provide clients $A(\u2208V)$ and $B(\u2208V)$ with maximally entangled state $|\Phi dkl\u27e9AB$ (or private state $\gamma \u0302dklAB$) with probability $pkl$ and error $\epsilon >0$ by using quantum channels $ne$$nqe$ times on average at most (i.e., $l\xafe\u2264nqe$), and $C$ represents a network quantum capacity $Q$ (private capacity $P$) per time. Notice that any protocol with $lr\xafe\u2264nqe$ belongs to the set $\Lambda (n,{qe}e\u2208E,\epsilon )$ of protocols.

We can capture alternative definitions of capacity by putting additional constraints on ${qe}e\u2208E$. 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: $\u2211e\u2208Eqe=1$. Notice that any protocol with $l\xafEf\xafe\u2264nqe$ belongs to the set $\Lambda (n,{qe}e\u2208E,\epsilon )$ of protocols with $\u2211e\u2208Eqe=1$.

Second, we can consider the usage of the whole network as the unit resource. This network use metric can be captured by setting $qe=1$ for all $e\u2208E$. This quantity is equivalent to the one introduced by Pirandola^{46} for multipath routing protocols with a so-called flooding strategy.

Finally, we note that Pirandola also introduced the single path per network use capacity.^{46} This quantity is not known to correspond to a restriction on ${qe}e\u2208E$. For additional details on this quantity we refer to the original Reference 46 and to Reference 67 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 $|\Phi d\u27e9AB$ is a special example of the private states $\gamma d\u0302AB$.

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\u2192B}$. The first useful upper bound on the private/quantum capacity is derived by Takeoka *et al.*,^{47} with the use of the squashed entanglement $Esq$, which applies to any point-to-point quantum channel $nA\u2192B$. 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.*,^{33} with the use of the relative entropy $ER$ of entanglement, to obtain tighter bounds for teleportation simulable channels, able to close the gap with previously-derived lower bounds.^{218} 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 $nA\u2192B$ is a single-mode lossy bosonic channel with a transmittance *η* ($0\u2264\eta \u22641$), $ER(nA\u2192B)=\u2212log2(1\u2212\eta )$ while we only know $Esq(nA\u2192B)\u2264log2[(1+\eta )/(1\u2212\eta )]$. 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,^{49} with the use of the max relative entropy $Emax$ 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.*^{68}

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,^{45,46} Azuma *et al.*,^{35} and Rigovacca *et al.*^{48} 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 $VA(\u2282V)$ and a party having the full control over $V\u2216VA$. In particular, Pirandola derives^{45,48} an upper bound (1) with $E=ER$ for any quantum network composed of teleportation simulable channels, with the use of the relative entropy $ER$ of entanglement, by following the PLOB bound,^{33} while Azuma *et al.* provide^{35} an upper bound (1) with $E=Esq$ for any quantum network composed of arbitrary quantum channels, with the use of the squashed entanglement $Esq$, by following the TGW bound.^{47} To make Pirandola's relative entropy bound applicable to arbitrary quantum networks like the upper bound of Azuma *et al.*, Rigovacca *et al.* present^{48} 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.^{49} This bound corresponds to the case where $E(ne)$ in Eq. (1) is assumed to be

where $ne\u2208S$ indicates that quantum channel $ne$ is teleportation simulable.

Pirandola gives^{45,46} an achievable protocol that works over a quantum network composed of teleportation simulable channels, while Azuma *et al.* present^{66} 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 $ne$. In particular, they assume that there is an entanglement generation protocol (including entanglement purification) over quantum channel $ne$ which provides a state $\rho \u0302e$ *δ*-close to $leR\delta e$ copies of a qubit maximally entangled state $|\Phi 2\u27e9AB$, called a Bell pair, by using quantum channel $(ne)\u2297le$, that is, $||\rho \u0302e\u2212(|\Phi 2\u27e9\u27e8\Phi 2|e)\u2297leR\delta e||1\u2264\delta $ with $\delta \u22650$ (for large *l ^{e}*). If this entanglement generation protocol runs over every channel $ne$, the network can share the state $\u2297e\u2208E\rho \u0302e$ with

If we regard each Bell state $|\Phi 2\u27e9\u27e8\Phi 2|e$ of the Bell-pair network $\u2297e\u2208E(|\Phi 2\u27e9\u27e8\Phi 2|e)\u2297leR\delta e$ as an undirected edge $e\u2033$ with the same two ends of *e*, we can make a multigraph $G\u2033=(V,E\u2033)$ with the set of $E\u2033$ of such undirected edges $e\u2033$. In this multigraph $G\u2033$, 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^{69,70} 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 Sec. II C 1 for this theorem). Thus, $M=minVA\u2211e\u2208\u2202(VA)leR\delta e$ for the multigraph $G\u2033$. 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 $minVA\u2211e\u2208\u2202(VA)leR\delta e$ Bell pairs with error $|E|\delta $. For fixed $qe=le/l(>0)$ (or $qe=le/lE(>0)$ with $lE=\u2211e\u2208Ele$), if the assumed entanglement generation protocol over every channel $ne$ can achieve quantum capacity, that is, $\delta \u21920$ and $R\delta \u2192Q(ne)$, by regarding *l* (*l ^{E}*) as

*n*and then taking $n\u2192\u221e$, the aggregated quantum repeater protocol gives lower bounds on the network quantum capacities of protocols $\Lambda (n,{qe}e\u2208E,\epsilon )$ as

If a quantum network is composed only of quantum channels $ne$ that satisfy Eq. (9), as Pirandola has considered,^{45,46} lower bounds in Eqs. (13) and (14) equal to upper bounds in Eqs. (5) and (6) with $E=ER$, 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.^{66}

#### 2. 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., $\tau \u0302dC:=\u2297i|\Phi d(i)\u27e9\u27e8\Phi d(i)|AiBi$ or $\tau \u0302dC:=\u2297i\gamma \u0302d(i)AiBi$, where $Ai,Bi\u2208V,\u2009C=\u2297iAiBi$ and $d:=(d1,d2,\u2026)$. The upper bounds for this problem are given^{45,71} by Pirandola for any quantum network composed of teleportation simulable channels, with the use of the relative entropy $ER$ of entanglement, while they are given by Bäuml and Azuma^{72} and Azuma *et al.*^{35} for any quantum network composed of arbitrary quantum channels, with the use of the squashed entanglement $Esq$. The works of Pirandola and Azuma *et al.* still use bipartite entanglement measures, while Bäuml *et al.* enable us to use multipartite squashed entanglement.^{73,74} The upper bounds are described as follows if we use the bipartite entanglement measure $E$ holding the bound (1): any protocol which provides every pair $AiBi$ of clients with maximally entangled state $|\Phi dkl(i)\u27e9AiBi$ (or private state $\gamma \u0302dkl(i)AiBi$) with probability $pkl$ and error $\epsilon >0$ (i.e., $||\rho \u0302klC\u2212\tau \u0302dklC||1\u2264\epsilon $), by using a quantum network associated with graph $G=(V,E)$, follows

for any $V\u2032\u2282V$, where $\u27e8xk\u27e9k:=\u2211kpkxk,\u2009\u2202(V\u2032)(\u2282E)$ is the set of edges which connect a node in $V\u2032$ and a node in $V\u2216V\u2032$, and $I(VA)$ is the set of indices *i* whose corresponding pair $AiBi$ satisfies $Ai\u2208V$ and $Bi\u2208V\u2216V\u2032$, or $Ai\u2208V\u2216V\u2032$ and $Bi\u2208V\u2032$.

Capacities for multiple-pair communication could have more variety than those for two-party communication because there are many pairs of clients.^{67} A possible definition is

for given ${qe}e\u2208E$ with $qe\u22650$, where $n\u22650,\u2009\Lambda (n,{qe}e\u2208E,\epsilon )$ is a set of protocols which provide clients $Ai(\u2208V)$ and $Bi(\u2208V)$ with maximally entangled state $|\Phi dkl(i)\u27e9AiBi$ (or private state $\gamma \u0302dkl(i)AB$) with probability $pkl$ and error $\epsilon >0$ (i.e., $||\rho \u0302klC\u2212\tau \u0302dklC||1\u2264\epsilon $), by using quantum channels $ne$$nqe$ times on average at most (i.e., $l\xafe\u2264nqe$), and $Cworst$ represents a worst-case network quantum capacity $Qworst$ (private capacity $Pworst$) per time. Notice that any protocol with $lr\xafe\u2264nqe$ belongs to the set $\Lambda (n,{qe}e\u2208E,\epsilon )$ of protocols. The use of objective function $mini{\u27e8\u2009log2dkl(i)\u27e9kl/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 $\u2211e\u2208Eqe=1$ on the set ${qe}e\u2208E$, the capacity per time is transformed into one per average total number of channel uses,

Notice that any protocol with $l\xafEf\xafe\u2264nqe$ belongs to the set $\Lambda (n,{qe}e\u2208E,\epsilon )$ of protocols with $\u2211e\u2208Eqe=1$. Another possible definition is

for given ${qe}e\u2208E$ with $qe\u22650$ and given ${si}i$ with $si\u22650$, where $Cweight$ represents a weighted multipair network quantum capacity $Qweight$ (private capacity $Pweight$) per time. The use of objective function $\u2211isi\u27e8\u2009log2dkl(i)\u27e9kl/n$ means that a protocol is evaluated with the average of rates for client pairs $AiBi$. By putting an additional constraint $\u2211e\u2208Eqe=1$ on the set ${qe}e\u2208E$, another possible definition is

where $Cweight$ represents a weighted multipair network quantum capacity $Qweight$ (private capacity $Pweight$) per the average total number of channel uses. As an important special case of $Cweight$, we may also consider total throughput $Ctotal$ defined as

and

Equation (15) gives upper bounds on the capacities. Let $Ri:=\u27e8\u2009log2dkl(i)\u27e9kl/n$. If we multiply Eq. (15) by $n\u22121$ with $n\u22650$ and take the limit of $n\u2192\u221e$ and $\epsilon \u21920$, it gives upper bounds on *R _{i}* for every

*i*, on $Ri+Rj$ for every $i\u2260j$, on $Ri+Rj+Rk$ for every different

*i*,

*j*,

*k*,…, and on $\u2211iRi$, by changing the choice of $V\u2032$. For given ${qe}e\u2208E$, if we minimize every upper bound over possible choices of $V\u2032$ for it, the minimized upper bounds give a minimal polytope in the Cartesian coordinates $R:=(R1,R2,\u2026)$ which includes all possibly achievable points of vector $R$. Since objective functions in the definitions of $Cworst(G,{ne}e\u2208E,{qe}e\u2208E)$ and $Cweight(G,{ne}e\u2208E,{qe}e\u2208E;{si}i)$ are concave, maximization of the objective functions over vector $R$ in the minimal polytope is convex optimization, which provides upper bounds on $Cworst(G,{ne}e\u2208E,{qe}e\u2208E)$ and $Cweight(G,{ne}e\u2208E,{qe}e\u2208E;{si}i)$. If we further change ${qe}e\u2208E$ to obtain upper bounds on $Cworst(G,{ne}e\u2208E)$ and $Cweight(G,{ne}e\u2208E;{si}i)$, 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

^{67}a computationally efficient method to obtain upper bounds, which are explained in Sec. 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.^{66} In particular, by running entanglement generation over every quantum channel $ne$, the network shares a quantum state $\u2297e\u2208E\rho \u0302e$ *ϵ*-close to $\u2297e\u2208E(|\Phi 2\u27e9\u27e8\Phi 2|e)\u2297leR\delta e$ associated with an undirected multigraph $G\u2033=(V,E\u2033)$, where each undirected edge $e\u2033\u2208E\u2033$ with the same two ends of *e* corresponds to a Bell pair $|\Phi 2\u27e9e$ between the ends of $e\u2033$. Then, if we find a path between *A _{i}* and

*B*in the undirected multigraph $G\u2033$, we can provide a state close to a Bell pair $|\Phi 2\u27e9AiBi$ by performing entanglement swapping along the $AiBi$-path. Thus, depending on an objective function in a definition of a capacity, by finding out edge-disjoint paths between pairs $AiBi$ in the undirected mutligraph $G\u2033$ properly, we can give a state close to $\tau \u0302dC:=\u2297i|\Phi d(i)\u27e9\u27e8\Phi d(i)|AiBi$. Along this recipe, Bäuml

_{i}*et al.*provide

^{67}a computationally efficient way to obtain lower bounds on the capacities, which are explained in Sec. II D.

#### 3. For multipartite communication

The upper bound^{72} of Bäuml *et al.* with the use of multipartite squashed entanglement^{73,74} works for broader classes of networks by generalizing the idea of Seshadreesan *et al.*^{75} 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 $ne$ is a directed hyperedge which has a single tail $t(e)\u2208V$ and multiple heads $h(e)(\u2282V$). 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 $\tau \u0302dC:=\u2297i|GHZd(i)\u27e9\u27e8GHZd(i)|Ci$ or multipartite private states $\tau \u0302dC:=\u2297i\gamma \u0302d(i)Ci$, where $Ci\u2282V$. Then, the upper bound is described as follows with the use of multipartite squashed entanglement $EsqP$: any protocol which provides every set *C _{i}* of clients with GHZ state $|GHZdkl(i)\u27e9Ci$ (or private state $\gamma \u0302dkl(i)Ci$) with probability $pkl$ and error $\epsilon >0$ by using a quantum broadcast network associated with a directed hypergraph $G=(V,E)$, follows

for any partition $P=P1:P2:\cdots :Pk$ of the set *V*, where $\u27e8xk\u27e9k:=\u2211kpkxk,\u2009EtriP(\u2282E)$ is the set of hyperedges whose tail *t*(*e*) and heads *h*(*e*) all belong to one set of *P*, $nCi|P=0$ for the case of $|{l\u2208{1,\u2026,k}:Pl\u2229Ci\u2260\u2205}|<2$ and $nCi|P=|{l\u2208{1,2,\u2026,k}:Pl\u2229Ci\u2260\u2205}|$ for the case of $|{l\u2208{1,2,\u2026,k}:Pl\u2229Ci\u2260\u2205}|\u22652$. The multipartite squashed entanglement of broadcast channel $ne$ is defined by $EsqP(ne):=max|\psi \u27e9t(e)EsqP1\u2229(t(e)h(e)):\cdots :Pk\u2229(t(e)h(e))(ne(|\psi \u27e9\u27e8\psi |t(e)))$, where if $Pj\u2229(t(e)h(e))=\u2205$, we strip $Pj\u2229(t(e)h(e))$ from the partition of $EsqP1\u2229(t(e)h(e)):\cdots :Pk\u2229(t(e)h(e))$ in the right-hand side, and $|\psi \u27e9t(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 $\u27e8\u2009log2dkl(i)\u27e9kl$ corresponding to the size of $|\Phi d(i)\u27e9AiBi$ or $\gamma \u0302d(i)AiBi$ in the definitions of Eqs. (3), (4), and (16)–(19) as one associated with $|GHZd(i)\u27e9Ci$ or $\gamma \u0302d(i)Ci$ in the current scenario. 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 Sec. II B 2). Lower bounds on the capacities can be obtained^{72} by borrowing the idea of the aggregated quantum repeaters,^{66} similarly to the scenarios for multiple-pair communication. In this case, by running entanglement generation over every quantum broadcast channel $ne$, the network shares a quantum state $\u2297e\u2208E\rho \u0302e$ *ϵ*-close to $\u2297e\u2208E(|GHZ2\u27e9\u27e8GHZ2|e)\u2297leR\delta e$ associated with an undirected multihypergraph $G\u2033=(V,E\u2033)$, where each undirected hyperedge $e\u2033\u2208E\u2033$ with the same ends of *e* corresponds to a GHZ state $|GHZ2\u27e9e$ among ends of $e\u2033$. Then, if we find a Steiner tree spanning clients *C _{i}*, which is an acyclic subhypergraph connecting all vertices of

*C*, in the undirected multihypergraph $G\u2033$, we can provide a state close to a qubit GHZ state $|GHZ2\u27e9Ci$ by performing generalized entanglement swapping

_{i}^{76}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*in the undirected multihypergraph $G\u2033$, we can give a state close to $\tau \u0302dC:=\u2297i|GHZd(i)\u27e9\u27e8GHZd(i)|Ci$. 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 graph theory.

_{i}^{77}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

^{67}a computationally efficient method to obtain upper and lower bounds, which are explained in Sec. II D.

Recently, an alternative bound on the conference key and GHZ rates is provided by Das *et al.*^{78} In this work, the entire network is 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 optimization theory

A number of results provide upper and lower bounds on the rates at which entangled resources for the tasks described in Sec. II B can be distributed in a quantum network given by an arbitrary graph.^{35,45,46,48,66,67,71,72} 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 Sec. II D). In particular, we review concepts such as network flows, the max-flow min-cut theorem, multicommodity flows, Steiner trees, and the complexity of the various optimization 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 Sec. II D.

In the following, we consider an undirected, weighted graph $G\u2032=(V,E\u2032)$, i.e., a graph consisting of vertices $v\u2208V$, undirected edges ${vw}\u2208E\u2032$ where each edge is assigned a weight $c{vw}\u22650$. 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 Sec. II D.

#### 1. 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}\u2208E\u2032$, we can define edge-flows $fvw\u22650$ and $fwv\u22650$ such that

The quantities *f _{vw}* and

*f*could be interpreted as the amount of the abstract commodity flowing from vertex

_{wv}*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\u2208V$, which we call the source and target, respectively. We can define a flow from *s* to *t* as

while requiring that for every $w\u2208V$ such that $w\u2260s$ and $w\u2260t$ it holds

Equation (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 maximize the flow $fs\u2192t$ 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,

where the maximization is over $2|E\u2032|$ edge flows $fvw\u22650$ and $fwv\u22650$. The LP given by Eq. (27) further has $|E\u2032|$ inequality constraints and $|V|\u22122$ equality constraints. Using $|E\u2032|$ slack variables, the $|E\u2032|$ inequality constraints can be transformed into equality constraints, resulting in an LP in a standard form with $N=3|E\u2032|$ non-negative variables and $M=|E\u2032|+|V|\u22122$ equality constraints. Using interior point methods,^{79} such an LP can be computed using $O(NL)$ iterations and $O(N3L)$ arithmetic operations, where *L* scales as $O(MN+M+N)$.^{80}

Next, let us introduce the concept of cuts. Given a subset $W\u2282V$ 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\u2216W$. In the case where $s\u2208W$ and $t\u2208V\u2216W$, we call Eq. (28) an *st*-cut and use the notation $\u2202(Ws: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,^{81,82} which states that for all $s,t\u2208V$ it holds

Let us note that so far we have only made the assumption that all edge flows are non-negative real numbers. Depending on the network scenario considered, however, it might be necessary to require the values *f _{vw}* and

*f*to be integers for all ${vw}\u2208E\u2032$. 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 {

_{wv}*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

^{69}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 optimization problem given by Eq. (27) into an integer linear program. Integer linear programing has been shown to be an NP-complete problem.

^{83}

#### 2. Multicommodity flows

The theory of network flows can be generalized to scenarios with *k* source-target pairs $(s1,t1),\u2026,(sk,tk)$, and *k* different commodities flowing through a network concurrently. To describe such a scenario, for every edge, we define 2*k* edge flows $fvw(i)\u22650$ and $fwv(i)\u22650$ for $i=1,\u2026,k$. We then require that

We can now define the flow of commodity *i* from *s _{i}* to

*t*as

_{i}Further we require that every commodity *i* is conserved in every vertex except its source *s _{i}* and target

*t*, i.e., for all $w\u2208V$ such that $w\u2260si$ and $w\u2260ti$, it holds

_{i}When maximizing 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 maximize the sum over the *k* flows, resulting in the maximum total concurrent flow

where the maximization is under the constraints that Eq. (32) holds for all edges ${vw}\u2208E\u2032$ and Eq. (34) holds for all commodities $i\u2208{1,\u2026,k}$ and for all vertices $w\u2208V$ such that $w\u2260si$ and $w\u2260ti$. Equation (35) is a linear program, which in a standard form has $N=(2k+1)|E\u2032|$ variables and $M=|E\u2032|+k(|V|\u22122)$ 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 maximized by a flow instance where the $fsi\u2192ti$ vary greatly. The second figure of merit we consider here is fairer in the sense that it maximizes the least flow which can be achieved for any of the commodities concurrently. Namely, we have the worst case concurrent flow

where again the maximization is under the constraints that Eq. (32) holds for all edges ${vw}\u2208E\u2032$ and Eq. (34) holds for all commodities $i\u2208{1,\u2026,k}$ and for all vertices $w\u2208V$ such that $w\u2260si$ and $w\u2260ti$. Adding a slack variable $fworst\u22650$, Eq. (36) can be reformulated as the maximization of $fworst$ such that, in addition to the above constraints, for all $i\u2208{1,\u2026,k}$ it holds $fworst\u2264fsi\u2192ti$, which is a linear program. This in a standard form has $N=(2k+1)|E\u2032|+1+k$ non-negative variables and $M=|E\u2032|+k(|V|\u22122)+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 maximization given by Eq. (27) in the case of a single commodity, i.e., *k* = 1. An obvious question arising now is whether the max-flow 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 generalize the minimum *st*-cut to multiple source–target pairs is to consider multicuts. For *k* commodities, $(s1:t1),\u2026,(sk:tk)$-multicut is a set of edges, the removal of which disconnects *s _{i}* from

*t*for all $i\u2208{1,\u2026,k}$. As in the case of cuts, the weight of a multicut is defined as the sum of weights of the edges it contains. Minimizing over all multicuts, we obtain the minimum multicut $cminmulticut$. In the case where

_{i}*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

^{84}and can be computed in polynomial time. However, in general, the problem of finding the minimum multicut is known to be NP-hard.

^{85}

A second way in which the minimum *st*-cut can be generalized to multiple source–target pairs is the min-cut ratio, which is defined as

where the minimization is over arbitrary subsets $W\u2282V$ and $d(\u2202(W))$ is the number of source–target pairs separated by $\u2202(W)$. Computation of the min-cut ratio has been shown to be NP-hard in general, as well.^{86}

Whereas it is clear from a complexity point-of-view that, for general *k*, there cannot be an equality relation between $cminmulticut$ or $cmincut\u2009ratio$ to a multicommodity flow instance that can be computed in polynomial time, there are relations up to a factor of $O(log\u2009k)$. Namely, it has been shown by Garg *et al.*^{87} that

where $g1(k)$ is of $O(log\u2009k)$. Similarly, it has shown by Aumann and Rabani^{86} that

where $g2(k)$ is of $O(log\u2009k)$.

#### 3. Steiner trees

When considering scenarios involving groups consisting of more than two network users rather than pairs of users, we need to generalize 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\u2033=(V,E\u2033)$ and a set $S\u2282V$ of vertices, a Steiner tree spanning *S*, or in short *S*-tree, is a connected, acyclic subgraph of $G\u2033$ connecting all $si\u2208S$. By definition, in the *S*-tree, any pair of vertices $si,sj\u2208S$ with $i\u2260j$ 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\u2033$. This problem, which is known as Steiner tree packing, is NP-complete in general.^{88} 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\u2033$: As in an *S*-tree any pair $si,sj\u2208S$ with $i\u2260j$ of vertices has to be connected by a path, the number $tS(G\u2033)$ of edge-disjoint *S*-trees in $G\u2033$ cannot exceed the number of edge-disjoint paths between any pair $si,sj\u2208S$ with $i\u2260j$ in $G\u2033$. Minimized over all $si,sj\u2208S$ with $i\u2260j$, this is known as the *S*-connectivity $\lambda S(G\u2033)$ of the graph. Application of Menger's theorem to each pair $si,sj\u2208S$ with $i\u2260j$ individually shows that $\lambda S(G\u2033)$ is equal to the minimum *S*-cut, which is defined as

Further, there have been a number of results^{89–91} providing lower bounds on the number of edge disjoint *S*-trees $tS(G\u2033)$, which are of the form

### D. Efficient bounds for communication tasks

As was shown by Bäuml *et al.*,^{67} a combination of the results presented in Secs. II B and II C 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 multicommodity flow optimizations 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 Sec. II B, we define an undirected graph $G\u2032=(V,E\u2032)$ as follows: if there is only a single directed edge connecting two vertices $v,w\u2208V$, we regard the edge as an undirected edge of the graph $G\u2032$ with the same two ends; if there exist two edges $vw\u2208E$ and $wv\u2208E$ for two vertices $v,w\u2208V$, we only add one undirected edge ${v,w}={w,v}$ to $E\u2032$. The reason we switch from the directed graphs to undirected ones here is that, whereas the original quantum channels $ne$ for $e\u2208E$ are directed, in protocols such as the aggregated repeater one presented in Sec. 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}\u2208E\u2032$ of $G\u2032$ can be equipped with two weights $cE({vw})$ and $cQ({vw})$ to be used in the flow maximisations for the upper and lower bounds, respectively. For protocols in $\Lambda (n,{qe}e\u2208E,\epsilon )$, we define

where we set $E(nvw)=Q(nvw)=0$ whenever $vw\u2209E$.

#### 1. 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\u2032$ with weights $cE({vw})$ for all ${vw}\u2208E\u2032$. By the max-flow min-cut theorem, Eq. (31), it can be expressed as a flow maximization of the form Eq. (27), which can be computed by linear programming. As the constraints in the maximization 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 maximization in the undirected graph $G\u2032$ with weights $cQ({vw})$, and hence the optimization in Eq. (14) becomes a linear program. In summary, we have^{67}

where we have defined $fQmaxA\u2192B(G\u2032)$ and $fEmaxA\u2192B(G\u2032)$ as maximum flow from *A* to *B* in $G\u2032$ with edge weights given by Eqs. (43) and (42), respectively. Further, the maximization is over $qe\u22650$ such that $\u2211e\u2208Eqe=1$ on both sides. Both the upper and lower bounds are linear programs of the form (27), with additional optimization 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.^{67} Roughly speaking, assuming that the $fQmaxA\u2192B(G\u2032)$ is maximized by noninteger edge flows *f _{vw}*, we can always find a large enough number

*N*such that

*Nf*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 $\u230aNfvw\u230b$ Bell pairs for all edges ${vw}\u2208E\u2032$, which can then be connected by entanglement swapping along paths linking Alice and Bob.

_{vw}#### 2. Multipair communication

Let us consider the case of multiple user pairs $(A1,B1),\u2026,(Ak,Bk)$ wishing to obtain bipartite resources such as Bell pairs or private states between them concurrently. In Sec. II B, it was shown that there are a number of different ways to define the corresponding capacities. One is to maximize the sum of the rates achievable by the user pairs. The corresponding quantum and private multipair network capacities $Qtotal$ and $Ptotal$, respectively, are defined in Eq. (19) with *s _{i}* = 1 for all

*i*. A second way is to maximize the worst case rate achievable between all user pairs, with corresponding network capacities $Qworst$ and $Pworst$, defined in Eq. (17).

It follows as a special case from the results of Bäuml and Azuma^{72} that the total private multipair capacity $Ptotal$ is upper bounded by the minimum multicut $cEsqminmulticut(G\u2032)$. Similarly, it has been shown^{67,72} that the worst-case private multipair capacity $Ptotal$ is upper bounded by the minimum cut ratio $cEmincut\u2009ratio(G\u2032)$ with respect to weights defined by Eq. (42) for all entanglement measures $E$ holding the bound (1). Applying the generalizations of the max-flow min-cut theorem to multicommodity flows, given by Eqs. (38) and (39), respectively, we can upper bound both private multipair capacities in terms of multicommodity flow optimizations.

Using a generalization of the aggregated repeater protocol,^{66} lower bounds on the multipair quantum capacities $Qtotal$ and $Qworst$ can be achieved by the corresponding multicommodity flow maximizations $fQmaxtotal(G\u2032)$ and $fQmaxworst(G\u2032)$ with respect to weights defined by Eq. (43).^{67} Again we note that the integer constraint in the flow optimization can be dropped in the asymptotic case. In summary, it holds for the total multipair network capacities

where the maximization is over all $qe\u22650$ with $\u2211eqe=1$. Both upper and lower bounds are linear programs. Similarly, for the worst-case multipair 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 discussed in Sec. II C 2, the gaps $g1(k)$ and $g2(k)$, both of which are of order $O(log\u2009k)$, leave the possibility for more efficient protocols, e.g., quantum network coding.

Let us also note that the multicommodity flow based approach to obtain achievable rates in multipair 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.*^{92} overcame this problem by adding an additional constraint on the lengths of paths to the optimization of the total multicommodity flow. They provide a linear program scaling polynomially with the graph parameters, which provides a maximal solution if it exists.

#### 3. Multipartite communication

Finally, we consider the case where a set of users $A={A1,\u2026,Ak}$ wishes to establish a *k*-partite GHZ or private state. It is straightforward to generalize 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 $QA(G,{ne}e\u2208E)$ and $PA(G,{ne}e\u2208E)$, respectively. From the results of Bäuml and Azuma,^{72} it follows that the multipartite private network capacity is upper bounded by the minimum $A$-cut in $G\u2032$ with weights $cEsq({vw})$,

By application of the max-flow min-cut theorem (for every pair $Ai,Aj\u2208A$ with $i\u2260j$), as well as introduction of a slack variable, the right-hand side can be expressed as a linear program.^{67} Namely, we maximize $f\u22650$ such that for all pairs $Ai,Aj\u2208A,\u2009i\u2260j,\u2009f\u2264fAi\u2192Aj$, 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 $QA(G,{ne}e\u2208E)$, we have to generalize 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\u2033$, remains the same. The second part, however, has to be generalized 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 generalized entanglement swapping protocol,^{76} resulting in GHZ states among $A$.^{67,72,93} Thus the relevant figure of merit is the number of edge disjoint $A$-trees in $G\u2033$.

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 optimization, which in the asymptotic limit can be further relaxed to a flow optimization problem of the form $minAi,Aj\u2208A,i\u2260jfQmaxAi\u2192Aj$, which is a linear program. In summary, we have the following bounds:

where we have set $g3=12$. Both the upper and lower bounds are linear programs with $2(|A|2)|E\u2032|+|E|+1$ variables and $(|A|2)(|E\u2032|+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.*,^{94} Schoute *et al.*,^{95} and Chakraborty *et al.*^{96} We will now give a brief, and by no means exhaustive, overview of other quantum network protocols, which go beyond the classical routing paradigm.

#### 1. Quantum network coding

The lower bounds presented in Sec. II D 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.^{66} 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.^{97}

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 multipair 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.^{98} 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.^{99}

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.^{100,101} 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.*^{100,101} to distribute Bell pairs concurrently between both user pairs. The protocol mainly consists of creation of the state $|+\u27e9=(|0\u27e9+|1\u27e9)/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.^{102}

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 multipair scenario.^{103} Note that there is a trade-off between the increased need of quantum memory in time-sharing and the need to perform more complex operations in quantum network coding.

#### 2. 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.*,^{104} 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 $|+\u27e9$-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 $|+\u27e9$-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 $|+\u27e9$-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 stabilizer codes to compensate for losses in the transmission channels, noisy gates, as well as errors in the preparation and measurement. In subsequent work^{105} the technique is combined with quantum network coding and generalized to qudit graph states.

The second task has been considered by Hahn *et al.*,^{106} 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^{107} 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.^{108} In related work by Dahlberg *et al.*,^{109,110} 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.

#### 3. Coherent control

The adaptive protocols described in Sec. 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 $n1$ and $n2$. By entering a state in quantum superposition, say $|+\u27e9c$ into the control register, a superposition of the orders of $n1$ and $n2$ is achieved.^{111}

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.^{112–114} In particular, it has been shown that two identical copies of a completely depolarizing 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.^{115}

Miguel-Ramiro *et al.*^{116} 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 Sec. II B.

## III. WAITING TIMES AND FIDELITY ESTIMATION OVER ABSTRACT QUANTUM NETWORKS

In this section, we review the 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.^{117} 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 Sec. III A with the modeling of a quantum network, which includes the mathematical abstraction of different components in a quantum network. In Sec. 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 Sec. III C: Markov chain methods, probability-tracking algorithms, and sampling with Monte Carlo methods. Finally, we review in Sec. 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.”

#### 1. 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,^{117} 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 $pg$. 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 $\rho \u2041$ that is produced is constant, i.e., it is independent of the number of attempts required to produce it. The state $\rho \u2041$ is a noisy Bell state which typically incorporates different sources of noise, photon loss, and detector inefficiency.

#### 2. Entanglement swapping

The photon transmission probability in fiber decays exponentially and thus fundamentally limits that distance over which elementary link generation can be performed.^{47} 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.^{118} Typically, entanglement swaps are probabilistic, with a fixed success probability $ps$ 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., $ps=1$. If entanglement swapping is implemented with optical components, the entanglement swapping becomes probabilistic, i.e., $ps<1$ and typically $ps\u22640.5$.^{119} There are also more sophisticated optical swapping schemes with a probability larger than one half.^{120–122} In some models, where the memory decoherence to the vacuum state is considered, the success probability can also be a variable.^{123}

#### 3. 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.^{124} In contrast with entanglement swapping, the success probability $pd$ depends on the entangled states that are distilled.^{124,125}

#### 4. Entangled state representation

Arguably, the simplest model of the fresh elementary link state is a Werner state,^{126} which characterizes the state with a single parameter *w*,

where $|\Phi 2\u27e9$ is the desired maximally entangled two-qubit state and $I\u0302/4$ is 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.^{28} A more general model is a probabilistic mixture of the four Bell states. This state is also always achieved^{56} by applying local Pauli gates randomly to an arbitrary state of a pair of qubits, and it is not less entangled than a Werner 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 $\rho \u0302$ of a state, its fidelity with a pure target state $|\varphi \u27e9$ is given by $\u27e8\varphi |\rho \u0302|\varphi \u27e9$. Throughout the section, the target state will be a Bell state.

#### 5. 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.^{127}

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.

#### 6. 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.

#### 7. Cutoff

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 cutoffs: if a link remains idle for too long, it is discarded. By discarding entanglement whose storage time exceeds some prespecified 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.^{128,129}

Note that the inclusion of cutoffs in entanglement distribution schemes complicates their analysis because of the additional effect of waiting time on the state quality.

#### 8. 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.^{130–137} In this section, unless explicitly mentioned, we follow the BDCZ scheme,^{130} 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 $2n$ ($2n+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 $2n$ 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 $2n$ segments (see Sec. III A) with only entanglement swapping, i.e., no distillation or cutoffs unless explicitly mentioned. For simplicity, we assume that the generation probability $pg$ is the same for each pair of adjacent nodes and the swapping probability $ps$ 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 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.

#### 1. 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/pg$ at the bottom level, $1/ps$ at the higher levels), yielding $Tn=1/psnpg$.^{130,131} 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 $Tn=(3/2)n/psnpg$.^{133,136,138–143} 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 Sec. III B 3.

The mean-only approach is a good approximation when $pg$ and $ps$ are much smaller than 1.^{138,139} 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 cutoffs. With this method, memory decoherence is either approximated by an inefficiency constant^{142} or studied only for the communication time.^{144} To provide a better estimation, one needs to consider the waiting time distribution and the statistics it results in, which we discuss below.

#### 2. 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 Secs. III B 3 and III B 4, 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\u2208{1,\u20092,3.\u2009..}$ and $p=pg$. 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 $pg\u226a1$, the geometric distribution can be approximated by an exponential distribution

which is a continuous distribution with $t\u22650$. Note that we have set the attempt duration Δ of entanglement generation to 1 (Sec. 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 $T0\u2032$ 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^{140}

After two elementary links are prepared, the repeater node performs an entanglement swap, which is a probabilistic operation with success probability $ps$. 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 $M0(k)$ are independent copies of *M*_{0}. Equation (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=ps$. 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 $\u27e8K\u27e9$ swap attempts until the first successful swap, and for each swap attempt, $\u27e8M0\u27e9$ 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 $pg<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 $qg=1\u2212pg$. The probability that one link is prepared *j* steps before the other is given by^{141,145}

Here we assume that $T0>T0\u2032$. Modeling the fidelity decrease as exponential-decaying function of the storage time, the fidelity of the earlier link decays by a factor $\Gamma =\u27e8\u2009exp\u2009(\u2212|T0\u2212T0\u2032|/tcoh)\u27e9$, where $tcoh$ 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^{141,145}

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 cutoffs^{128} or protocols where two elementary links have to be prepared sequentially.^{146}

Unfortunately, for higher-level nested protocols, i.e., $n\u22651$, there is no analytical expression known for the mean waiting time $\u27e8Tn\u27e9$ with $ps<1$ because *T _{i}* for

*i*>0 does not follow a geometric distribution, in contrast to

*T*

_{0}.

#### 3. 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 $\u27e8Tn\u22121\u27e9$ and then approximate the distribution of $Tn\u22121$ with a geometric distribution [Eq. (49)] parameterized by $p=1/\u27e8Tn\u22121\u27e9$. Under this assumption, the mean waiting time $\u27e8Tn\u27e9$ can be computed by a derivation analogous to the one leading to Eq. (55) in Sec. III B 2 and is given by

with $pn\u22121=1/\u27e8Tn\u22121\u27e9$. In the limit of $pg\u21920$ for the bottom level (*n* = 0) and $ps\u21920$ for higher levels (*n* > 0), the mean waiting time $\u27e8Tn\u22121\u27e9\u2192\u221e$ and thus $pn\u22121\u21920$. 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 $\u27e8T0\u27e9=1/pg$ yields

which is precisely the 3-over-2 approximation mentioned in Sec. III B 1.

The error introduced by the approximations Eqs. (60) and (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., $pg\u21920$ and $ps\u21920$. 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 $ps$ 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 $pg$ with $1/\u27e8Tn\u22121\u27e9$. 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 decoherence 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 semianalytical method for computing fidelity with more sophisticated memory decay models^{148,149} (see also Sec. 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.*^{150} introduced 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 cutoff (see Sec. III A); with a memory cutoff, 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.* showed 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.^{151} 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 prespecified time with high probability.

#### 4. 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 ($ps=1$) and without distillation.

First, let us compute the waiting time probability distribution in the case of deterministic swaps without distillation or cutoffs. 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. For a repeater scheme with *n* nesting levels, the total waiting time is the time until all $N=2n$ elementary links are prepared, 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 $TN,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(TN,k\u2264t)$, 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^{152}

where $Pr(T0\u2264t)=1\u2212(1\u2212pg)t$ since *T*_{0} is geometrically distributed with success probability $pg$. 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 $\u27e8TN,k\u27e9$ from Eq. (62) can be computed by solving a recurrence formula where $\u27e8TN,k\u27e9$ is expressed as a function of $\u27e8TN,k\u22121\u27e9$.^{153,154} For *k* =* N*, i.e., the waiting time that all elementary entanglement are prepared,^{153} the solution reads

For $pg\u226a1$, this expression can be approximated by^{145,147,152}

with

where *H*(*N*) is the *N*th harmonic number and $\gamma \u22480.57721$ is the Euler–Mascheroni constant. In separate work, Praxmeyer included finite memory time with cutoff into the calculation and obtained^{154}

where *τ* is the cutoff threshold and $qg=1\u2212pg$.

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 Refs. 128 and 152.

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.^{152,153}

These expressions presented here for the deterministic-swap 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.^{155} 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 elementary-link success probability $pg$.

The waiting time in the deterministic-swap case can be used as an approximation to the case where $ps$ is slightly lower than 1 and bounds from below the waiting time for general $ps$. The quality of the approximation is shown in Fig. 4.

#### 5. 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 $pd$ 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.^{125,130} 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^{130} and memories.^{144} 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 cutoff, i.e., protocols that discard the entanglement if the fidelity is lower than a certain threshold.^{156} Because the distillation success probability is a monotonic function of the fidelity of the input states, a lower bound on the fidelity by a cutoff also directly yields a lower bound on the success probability.

### C. Numerical tools for evaluating analytical expressions

Above, in Sec. 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, closed-form expressions is unfeasible. An example of such a case is a nested repeater chain with cutoffs and nondeterministic 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 Sec. III B 2, where the single-repeater 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 Sec. IV in that they are all directly applied to expressions that can in principle be evaluated analytically.

#### 1. 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,^{147} 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\xaf$) occurs with the entanglement swapping success probability $ps$ (entanglement swapping succeeded), whereas the transition to “no entanglement” (state 00) has probability $1\u2212ps$ (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\xaf$ in Fig. 5). Shchukin *et al.* used this equivalence to compute the average waiting time, as well as its variance, by solving a linear equation system that has the same size as the number of states in the Markov chain.

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 and Kok.^{156} 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,^{130} 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 Sec. 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.*,^{157,158} 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.^{157}

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 reduced. With this technique, Shchukin *et al.* gave examples for the BDCZ protocol with analytical expressions for up to four nodes, while numerically they reached 32 nodes.^{147} Vinay and Kok reduced the computational complexity of this approach using probability generating functions and complex analysis, but the scaling remains exponential.^{156} To process a larger number of nodes, Shchukin *et al.* proposed to use the average waiting time to replace the random variable for low-level subprotocols. This idea is similar to approximating the waiting time distribution at every nesting level of a repeater protocol with the bottom-level distribution (see Sec. III B 3).

Finally, in a recent development, Khatri^{159} introduced a method for describing network protocols based on quantum partially observable Markov decision processes.^{160} 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 the protocol uses to choose its next action. As an application of the method, Khatri found analytical solutions for optimizing a cutoff 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.

#### 2. 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 $2n$ segments and no entanglement distillation [depicted in Fig. 2(b) for *n* = 2]. In Sec. 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 $Tn\u22121$ on the previous level. Unfortunately, as discussed in Sec. III B, to our knowledge, this recursive expression of

*T*has not been analytically evaluated for

_{n}*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 $ttrunc$, the waiting time distribution becomes finite. The evaluation with numerical tools leads to an algorithm that tracks both the truncated probability distribution of

*T*and the associated fidelity.

_{n}^{161,162}

On the $2n$-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,^{161} a computationally faster approach is to separate the probability distributions of failed and successful swap attempts.^{162} 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.^{162} 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 $ttrunc$ and has been used to track over 1000 nodes for some parameter regimes.^{161}

A related approach to the probability-tracking algorithm explained above is taken by Kuzmin *et al.*^{148,149} 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 memory 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.

#### 3. 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 semianalytical computation is to sample the expressions on random variables for the waiting time and the delivered state using a Monte Carlo approach.^{161} 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 $2n$-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 cutoffs. 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 one-way quantum error correction) and third-generation repeaters (loss-tolerant entanglement generation).^{117,163}

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,^{164–167} entanglement swapping and (one-way) quantum error correction are no longer probabilistic. (We note that although entanglement generation is also 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,^{163,168–178} 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 second and third generation) and the propagation of physical loss errors into logical errors (for third 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 in-depth discussion about network simulation, we refer to the books.^{179–181}

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 (Sec. III) or of the optimal achievable performance (Sec. 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. Refs. 1 and 179.

#### 1. 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.

#### 2. Theoretical modeling

The second step is modeling. Best practices^{1,179} 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. In particular, it 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.

#### 3. 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,^{179} a paradigm for simulation that guarantees reliable timing behavior (see Fig. 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 multithreaded 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.

#### 4. 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 cannot 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 cannot simulate computations involving more than tens of qubits.^{182,183}

Fortunately, many scenarios of interest in quantum networks involve only Clifford gates and Pauli noise, yielding classical simulation tractable.^{184} 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.^{185–194}

The first two platforms are SimulaQron^{195} and the QUantum NETwork SIMulator (QuNetSim).^{196} 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^{197} 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.^{198} 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^{199} and several multipartite protocols;^{200,201} moreover, within the context of the quantum protocol zoo,^{202} SimulaQron has been used to implement other protocols such as a quantum cheque scheme,^{203} a protocol for leader election,^{204} universal blind quantum computation,^{205} and coin flipping.^{206}

The next platform is the Simulator for Quantum Networks and CHannels (SQUANCH).^{207} 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^{208} and to investigate a photonic lattice architecture for a universal quantum programmable gate array.^{209}

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)^{210,211} aims to capture complex quantum network behavior. In particular, one of the target use cases is understanding the 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.^{212} This implies that most of OMNeT++ 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)^{213} and the NETwork Simulator for QUantum Information with Discrete events (NetSquid).^{214} 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^{215,216} as well as for analysis of a nine-node quantum network with quantum memories modeled after single erbium ions in solids.^{213} 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.^{214} NetSquid has also been used to evaluate a quantum router architecture^{217} and to go beyond the analytical regime in the analysis of a quantum switch.^{157,214}

## 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 multiuser 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 cutoffs, 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 validating protocols, evaluating alternative network solutions, studying 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.

## ACKNOWLEDGMENTS

This work was supported by the QIA project (funded by European Union's Horizon 2020, Grant Agreement No. 820445) and by the Netherlands Organization for Scientific Research (NWO/OCW), as part of the Quantum Software Consortium program (Project No. 024.003.037/3368). K.A. is thankful for the support, in part, from PREST, JST JP-MJPR1861 and from CREST, JST JP-MJCR1671. S.B. acknowledges support from the Government of Spain (FIS2020-TRANQI and Severo Ochoa CEX2019-000910-S), Fundació Cellex, Fundació Mir-Puig, Generalitat de Catalunya (CERCA, AGAUR SGR 1381) as well as from the European Union's Horizon 2020 research and innovation programme, Grant Agreement No. 820466 (project CiViQ). B.L. was supported by the IDEA League student grant programme.

## DATA AVAILABILITY

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