Utility-Based Path Selection and Configuration in Quantum Networks via Layered Shortest Paths Leonardo Bacciottini, Subhransu Maji, Don Towsley, Gayane Vardoyan
arXiv:2609.16198v1 [quant-ph] 14 Sep 2026
Manning College of Information and Computer Sciences University of Massachusetts Amherst Email: {lbacciottini,smaji,dtowsley,gvardoyan}@umass.edu
Abstract—A path in a quantum network is a chain of repeaters that distributes entanglement between two users. Selecting a path requires balancing the rate and quality (e.g., fidelity) of the delivered entanglement, but these quantities, unlike standard routing metrics, compose non-additively. The problem is compounded by link-level configuration choices (e.g., distillation rounds or emitter brightness tuning), each trading rate against fidelity, so that a path’s performance depends jointly on its route and its per-link settings. We cast this joint path selection and configuration problem as a shortest-path computation on a layered graph whose layers track discretized end-to-end fidelity. A single run returns the full rate–fidelity Pareto frontier, from which the path maximizing any nondecreasing utility function of rate and fidelity can be selected. We prove that for certain utility functions (including the secret key rate of BB84), the method is a fully polynomial-time approximation scheme, returning a nearoptimal path within a specified tolerance. We further characterize exactly when cheaper scalarization-based routing suffices: it is optimal for utility functions with convex fidelity profiles, but can be arbitrarily suboptimal otherwise (e.g., for step-like, sigmoidal utilities), whereas the layered method remains reliable in all cases.
I. I NTRODUCTION Quantum networks distribute entanglement to support applications such as provably secure communication [1], [2], quantum-enhanced sensing [3]–[6], verifiable blind quantum computation [7], and distributed quantum computation [8] with up to exponential speedups compared to its conventional counterpart. Quantum network performance depends critically on how entanglement is routed through chains of lossy and noisy quantum channels. Unlike standard classical routing metrics, however, the two principal measures of an entanglement path—the rate and quality (e.g., fidelity) of the delivered end-to-end states—compose non-additively and generally trade off against one another. Optimizing either quantity alone can therefore lead to a poor application-level solution. A high-rate path may produce states whose fidelity is too low to be useful, while a high-fidelity path may distribute entanglement too slowly. Application objectives such as the secret key rate of a quantum key distribution (QKD) protocol depend jointly on both quantities and are better described by a quantum utility function [9]. Selecting a path according to such a utility allows the routing decision to reflect the requirements of the application being executed, rather than optimizing a single intermediate performance measure such as rate, fidelity, or hop count.
The routing decision is further complicated by configurable rate–fidelity tradeoffs at individual links. For example, the bright-state population in single-click entanglement generation protocols [10], or the number of entanglement distillation rounds can be adjusted to trade generation rate for state quality (see more discussion in § III). We refer to the available configuration choices on a link as its operating points. Consequently, the performance of an end-to-end connection depends jointly on the selected physical route and the operating point chosen on every traversed link. We study this joint path-selection and link-configuration problem for bipartite entanglement distribution, where the target is an ideal Bell state to be shared between two nodes in the network. Given a pair of nodes s and t, and a nondecreasing utility function U (R, F ) of end-to-end rate R and fidelity F , our goal is to select an s–t path and one operating point on each link on the path to maximize the resulting utility. Our formulation builds on the quantum network utility maximization framework [9], but focuses on path selection and explicitly incorporates discrete per-link configuration choices. For clarity of presentation, we develop our formulation under the assumption of Werner states [11] at the elementary link level. Werner states are commonly used to model worstcase Pauli noise in quantum networks [12], thus enabling studies of performance lower-bounds assuming Pauli channels. A Werner state can be described in terms of its parameter W , which is linearly related to its fidelity F (see § III). Our key observation is that both rate and fidelity become additive properties after a suitable reparameterization. We exploit this structure though a layered graph whose layers discretize the end-to-end fidelity, while its edge costs track the corresponding rate. Running Dijkstra’s algorithm [13] on this graph computes the highest-rate path for every node and fidelity layer. This construction yields the quantum layered shortest-path (QLSP) algorithm. In contrast to methods that collapse rate and fidelity into a single edge cost, QLSP retains the end-to-end fidelity of the state explicitly and can therefore recover Pareto-optimal solutions that are not supported by any weighted-sum scalarization of functions of rate and fidelity. Moreover, given a node s, a single run computes the discretized frontier for every network node t and can be reused for any nondecreasing utility function. We analyze the approximation introduced by fidelity discretization and provide polynomial-time approximation guarantees under
certain assumptions (see Theorem 1 in § IV). We also establish a sufficient condition under which a weighted-sum scalarization recovers the utility-optimal route and link configuration (Theorem 2, § IV). This result motivates a computationally cheaper baseline that is provably optimal whenever this condition holds. This condition holds for the asymptotic secret key rate of the BB84 QKD protocol [1], as well as a utility based on entanglement negativity [14]. However, nonconvex fidelity profiles, such as sigmoid functions or those induced by hard fidelity requirements, can have Pareto-optimal but unsupported maximizers that no weightedsum scalarization can recover. QLSP continues to retain such solutions because it approximates the full rate–fidelity frontier rather than only its supported convex-hull vertices. Our main contributions are as follows: • We introduce QLSP, a layered shortest-path algorithm that solves the joint path selection and configuration problem by computing the rate–fidelity Pareto frontier from a given node s to every node t in a single shortest-path computation. • We analyze the complexity and approximation quality of QLSP, and establish polynomial-time approximation guarantees under the assumptions stated in § IV. • We characterize when a linear scalarization is exact, and propose a cheaper algorithm that exploits this property. We evaluate QLSP on both synthetic and real-world network topologies, including heterogeneous architectures. The experiments show that in the considered cases there are Paretooptimal candidate solutions that no linear scalarization can find. We also observe that QLSP outperforms all other considered baselines, including those that optimize only the endto-end fidelity or rate. II. R ELATED W ORK We organize prior work based on treatment of the rate– fidelity tradeoff. A comprehensive taxonomy of entanglement routing appears in the recent survey of Abane et al. [15]. a) Rate- and fidelity-aware quantum routing: Early work on repeater-network path selection ranks paths using additive per-link costs derived from simulated link throughputs [16], an approach representative of the edge-based baselines considered in § IV. Subsequent methods optimize end-to-end rate or network-wide throughput, including single-path, multipath, and concurrent-flow formulations [17]–[24]. Fidelity is typically ignored or imposed as a hard constraint. For example, Q-PATH/Q-LEAP jointly selects a path and per-link purification rounds subject to an end-to-end fidelity threshold [25], while later work jointly optimizes routing and purification to maximize throughput under a fidelity requirement [26]–[28]. Related scheduling methods allocate resources to satisfy rate and fidelity requirements on fixed paths [29]. Such fidelityconstrained objectives are special cases of our framework, represented by utilities such as U (R, F ) = R 1[F ≥ Freq ], where R is the rate of entanglement generation, F is the fidelity to the desired entangled state, and Freq is a prespecified fidelity threshold imposed on each distributed state.
b) Utility-based routing: Our objective builds on the QNUM framework of Vardoyan et al. [9], which adapts classical network utility maximization [30]–[32] to entanglement distribution over fixed routes, using either centralized [9] or distributed [33] optimization. Most relevant, Kar and Mukhopadhyay [34] jointly optimize routes and continuous rate–fidelity allocations through a mixed-integer convex formulation. They obtain exact solutions for negativity-based utilities in a high-rate regime, and otherwise use a convex relaxation. Using a discrete set of rate–fidelity configurations, Zhang et al. [35] jointly optimize path selection, per-link configuration, and purification for multiple rate demands under a hard end-to-end fidelity constraint, solved with a greedy shortest-path heuristic refined by Bayesian optimization. Our setting considers a single pair of end nodes wishing to share entanglement, but allows arbitrary monotone utility functions and comes with approximation guarantees unlike theirs. c) Multi-criteria shortest paths in classical networking: Routing under two additive metrics is a classical problem in graph theory. In particular, finding a minimum-cost path subject to a delay constraint is NP-hard [36], admits fully polynomial-time approximation schemes (FPTASs) based on rounding-and-scaling dynamic programs [37], [38], and is often addressed in practice using Lagrangian relaxation, most notably the LARAC algorithm [39]. LARAC efficiently searches over the multiplier space, but it inherits the duality gap associated with scalarization; moreover, the number of supported solutions it may need to explore can be superpolynomial [40]. This also renders techniques that maintain the pareto-frontier explicitly more challenging (e.g., Coutinho et al. [41]), as the size of the frontier is not bounded a priori. QLSP adapts the rounding-and-scaling framework to quantum networks through the reparameterization of rate and fidelity which makes the quantum composition rules additive: rates compose harmonically (under certain assumptions, see § III), while fidelities compose multiplicatively. It further incorporates discrete perlink operating points, supports arbitrary monotone utilities, and provides an exact characterization of optimality (Theorem 2). III. M ODEL We consider first-generation quantum networks [42] which distribute entanglement to end nodes first by generating linklevel entanglement (LLE) between intermediate, physicallyadjacent network nodes, called quantum repeaters. Quantum repeaters then “fuse” these shorter-distance entangled states into end-to-end (e2e) entanglement via Bell-state measurements (BSMs). This process is also known as entanglement swapping. LLE generation (LLEG) is an inherently probabilistic process that degrades exponentially with link length in optical fiber: the success probability pgen ∝ exp{−L/Latt }, where L is link length and Latt is the fiber attenuation length, usually taken to be 22 km. Unless specified otherwise, we take all link lengths L to be in units of km. If the architecture allows it, distillation protocols [43], [44] can also be used at any stage to probabilistically to reduce the number of entangled states and improve quality.
Quantum noise model: Realistic quantum systems produce non-ideal quantum states that further degrade due to imperfect quantum storage, manipulation (gates), and measurement. Fidelity is a measure of closeness between quantum states (or gates) and their ideal realizations. We assume in this work that LLEG on link l produces a Werner state
(Rk , Wk ) denote the expected rate and Werner parameter after k rounds, so that k = 0 is the undistilled raw pair. Under ideal local operations and perfect storage of quantum states, the success probability for the next BBPSSW round is
1 − Wl I4 (1) 4 with probability pl , where Wl ∈√[0, 1] is the state’s Werner parameter, |Ψ− ⟩ = (|01⟩ − |10⟩)/ 2 is the desired Bell state, and I4 is the 4×4 identity matrix. The fidelity of ρl to |Ψ− ⟩ is then Fl = (3Wl + 1)/4, and for the state to be considered entangled it is necessary that Fl > 0.5 (equivalently, Wl > 1/3). While in reality LLEG can produce states that are not Werner, we assume their twirling [45] into the form (1) to simplify analysis. Namely, Werner states have the convenient property that entanglement swapping on a set of n Werner Qnstates ρl yields another Werner state with parameter W = l=1 Wl . Operating points: Physical mechanisms for entanglement generation commonly allow generation rate to be traded for state quality. This occurs, for example, when tuning the brightstate population in single-click LLEG [10], the pump power in spontaneous parametric down-conversion (SPDC) [46], or the number of rounds in an entanglement distillation protocol. Each available configuration therefore gives a different expected LLEG rate and Werner parameter. We call one such rate–quality pair (R, W ) an operating point. For single-click entanglement generation, the bright-state population α ∈ (0, 1) controls the rate–fidelity tradeoff. In the regime of small link transmissivity η, the fidelity and success probability are
and the output Werner parameter and expected rate obey [47]
ρl = Wl Ψ− Ψ− +
F (α) = 1 − α,
pgen (α) = 2η(L)α.
(2)
Assuming a heralding station placed exactly midway between two nodes separated by L km, the “link” is defined as the node-to-heralding station segment of length L/2 km, yielding η(L) = ce−L/(2Latt ) ,
(3)
where 0 < c < 1 collects coupling, conversion, detector, and other efficiencies not due to fiber attenuation. If LLEG is attempted at repetition rate Rrep , a setting α maps to the QLSP operating point R(α) = Rrep pgen (α),
W (α) = (4F (α) − 1)/3.
(4)
Thus, increasing α raises the generation rate while lowering the Werner parameter. Although α is continuous, selecting A admissible settings A α1 , . . . , αA produces the finite menu R(αa ), W (αa ) a=1 used by QLSP. Entanglement distillation produces the same algorithmic interface through a different physical tradeoff: it consumes multiple lower-quality pairs to produce fewer pairs of higher quality. Consider the BBPSSW protocol [43], which distills two identical Werner states into one. Suppose that LLEG supplies raw states with Werner parameter W0 at expected rate R0 , and allow at most n nested distillation rounds. Let
Psucc (Wk ) = (1 + Wk2 )/2,
Wk+1 =
2Wk (1 + 2Wk ) , 3(1 + Wk2 )
Rk+1 =
(5)
Psucc (Wk ) Rk . (6) 2
Choosing between zero and nn nested rounds therefore gives the finite menu (Rk , Wk ) k=0 . Unlike single-click tuning, the number of rounds is already a discrete protocol choice, so this entire menu is represented directly. The single-click scheme and BBPSSW distillation thus have different physical machinery but expose the same QLSP input: a finite collection of feasible (R, W ) pairs. Other physical or protocol controls can be incorporated in the same way. Sampling a continuous tradeoff more finely enriches the set of available configurations, at the cost of additional computation; this physical-menu discretization is distinct from the algorithmic approximation introduced in § IV. Utility functions: The QLSP algorithm can accommodate arbitrary quantum utility functions U (R, F )—where F is the fidelity of the average e2e Werner state distributed to two network nodes and R is the e2e generation rate—provided that U is monotonically increasing in its inputs. The utility function captures the complex ways in which users or quantum applications balance rate and fidelity. As an example, the asymptotic secret key rate (SKR) of the BB84 quantum QKD protocol, when carried out with Werner states of fidelity F distributed at rate R, reads o n 2(1 − F ) , (7) UBB84 (R, F ) = max 0, R 1 − 2h 3 where h(·) is the binary entropy function. Utility functions can combine any function of rate and an entanglement monotone, e.g., negativity, which for Werner states simplifies to n 1 o Uneg (R, F ) = max 0, R F − . (8) 2 Another commonly used one is the fidelity-threshold utility, Uth (R, F ) = R1[F ≥ Freq ]
(9)
which maximizes the generation rate subject to a fidelity requirement Freq . It may also be desirable to apply transformations to these utilities, such as a logarithm, to promote fairness among coexisting network flows [9]. Path selection under such nonlinear objectives is even more challenging when network links have their own configurable rate–fidelity tradeoffs, since optimizing each link according to a local objective does not in general produce an e2e optimum. These utilities can equivalently be expressed in terms of the Werner-state parameter W using the relation F = (3W +1)/4. In the remainder of the paper, we also use this parameterization because it simplifies the analysis.
Problem statement: Consider a quantum network represented by a graph G = (V, E), where vertices v ∈ V represent nodes (e.g., repeaters and end nodes), and undirected edges e ∈ E represent physical network links. Each edge e ∈ E may admit one or more operating points a ∈ Ae . Each operating point is associated with a rate-Werner parameter pair (Re,a , We,a ), where Re,a > 0 is the average rate at which LLEG produces a new entangled pair, and We,a ∈ [0, 1] is the corresponding Werner parameter. Given nodes s and t, producing an entangled state between them requires selecting a routing solution π = (p, a), where: 1) p = (e1 , e2 , . . . , en ) denotes a path from s to t, represented as a sequence of edges ek ∈ E; and 2) a = (a1 , a2 , . . . , an ) denotes the corresponding operating points, i.e., ak ∈ Aek is the operating point on edge ek . Each routing solution π = (p, a) induces an e2e rate Re2e (π), representing the average number of entangled states produced between s and t per unit time, and an e2e Werner parameter We2e (π), capturing their quality. Under the models adopted in this paper, they compose as Re2e (π) =
!−1
X
1
k
Rek ,ak
, We2e (π) =
Y
Wek ,ak . (10)
k
Re2e is the reciprocal of the expected e2e latency and is exact when entanglement swaps are performed sequentially along the path. A bottleneck model would instead use Re2e (π) = mink Rek ,ak . We use the first relation throughout, although our framework extends to any additive rate model. The second relation (We2e ) follows from Werner-state composition. Depolarizing gate noise can be included through additional multiplicative factors; memory decoherence can be included only when represented by a constant depolarizing term. Given nodes s and t, and a utility function U (R, W ) of rate R and Werner parameter W , our goal is to find a routing solution π = (p, a) that jointly selects and configures an s–t path to maximize its end-to-end utility: π ⋆ = arg max U (Re2e (π), We2e (π)). π
(11)
We note that it is possible to use time multiplexing or multipath routing to obtain solutions that use multiple paths or configurations to get a better objective than (11), but these techniques are out of the scope of this work. IV. M ETHODS In this section, we formally present the QLSP algorithm for approximately solving (11) and analyze its correctness, computational complexity, and approximation guarantee. We then introduce several baselines that we evaluate against QLSP in § V. Finally, in Theorem 2, we show that one of these baselines—the surrogate sweep algorithm—is optimal for quantum utility functions with convex fidelity profiles.
A. Quantum layered shortest-path (QLSP) formulation Because quantum-network utility functions may depend on both e2e rate and fidelity, directly optimizing either quantity alone does not generally yield an optimal solution. A path with the highest e2e rate, for instance, might yield zero utility if its e2e fidelity is too low (e.g., ≤ 1/2 for UNEG ). Our formulation (Algorithm 1) instead maintains a discretized Pareto frontier over the achievable rate–fidelity tradeoffs. Specifically, for each node and each discretized value of the e2e Werner parameter, we maintain the best e2e rate among paths terminating at that node. This construction can be viewed as introducing multiple layers for each node, with each layer corresponding to a discretized Werner-parameter value. We then show that computing shortest paths in this layered graph recovers the discretized Pareto frontier, which contains an approximately optimal solution, and we quantify the approximation error. a) Path composition rule: We assume that, compatibly with (10), extending a path by edge e operating at point a ∈ A changes the path’s effective rate R and Werner parameter W according to: 1 1 1 = + , W ′ = W We,a . (12) ′ R R Re,a It is convenient to reparameterize these quantities as: X ≡ 1/R, and Y ≡ − ln W . Under this change of variables, the update becomes additive: X ′ = X + xe,a , and Y ′ = Y + ye,a , where xe,a := 1/Re,a , and ye,a := − ln We,a . b) Layered graph construction: Choose a discretization width ∆ > 0 and let the Werner parameter bins be indexed by k = 0, 1, . . . , K. Bin k represents cumulative fidelity-loss Y approximately equal to k∆, corresponding to Werner parameter Wk ≈ e−k∆ . Given a Werner-parameter threshold Wmin , and K can be set as K = ⌈ln(1/Wmin )/∆⌉. For example, one may choose Wmin > 1/3 when optimizing negativity, since states below this threshold have zero negativity. The value Wmin is not intrinsic to the construction; it simply truncates the layered search space. We construct a layered graph G∆ = (V ∆ , E ∆ ) as follows. The layered node set is V ∆ = {(v, k) : v ∈ V, k ∈ {0, . . . , K}}. Thus, each physical node v is replicated across K + 1 fidelity layers (equivalently, Werner-parameter layers, since there is a one-to-one correspondence between them). For each physical edge e = (u, v) ∈ E, leach operating point a ∈ Ae , and each m k∆+ye,a ′ layer k, define k ′ = ≤ K, then we add a . If k ∆ ′ directed layered edge (u, k) → (v, k ) with nonnegative cost c (u, k), (v, k ′ ) = xe,a = 1/Re,a . If the physical graph is undirected, we add the reverse layered edge as well. We show an example of how the layered graph is constructed in Fig. 1. c) Interpretation: A path π = (p, a) in the layered graph corresponds to a physical path p together with a choice a of operating point on each traversed edge. The pathPcost in the layered graphP is the cumulative inverse-rate: X = e∈p 1/Re , and Y ≈ e∈p − ln We . Hence, if d(v, k) denotes the shortest-path distance from (s, 0) to (v, k) in the layered graph, then R(v, k) = 1/d(v, k), and W (v, k) ≈ e−k∆ .
Algorithm 1 The QLSP algorithm for rate–Werner utility
Physical Graph
𝑒!
s
u
𝑎" = 𝑅" , 𝑊" 𝑎! = 𝑅! , 𝑊! 𝑎# = 𝑅# , 𝑊#
𝑒"
v
𝑎" = 𝑅" , 𝑊" 𝑎! = 𝑅! , 𝑊!
Require: Graph G = (V, E), source s, destination t, fidelity bin width ∆, edge operating points {(Re,a , We,a )}a∈Ae . Ensure: Approximate path maximizing U (Re2e , We2e ) 1: Construct layered node set
t
Operating points
V ∆ = {(v, k) : v ∈ V, k = 0, . . . , K}.
Layered Graph s, 0
1/𝑅"
u, 0
t, 0
v, 0
…
…
v, 𝑘′"
…
…
u, 𝐾
v, 𝐾
t, 2 t, 3
…
u, 𝑘#
v, 𝑘′! 1/𝑅!
… …
Pareto frontier
u, 𝑘"
t, 1
1/𝑅"
…
… 1/𝑅#
…
u, 𝑘!
…
…
1/𝑅!
t, 𝐾
Fig. 1. Example of a physical and its associated layered graph. In the example, the edge e1 has three operating points a1 , a2 , a3 . In the layered graph, these correspond to three edges from (s, 0) to (u, k1 ), (u, k2 ), (u, k3 ), where ki = ⌈− ln Wi /∆⌉ with edge cost 1/Ri respectively. The edge e2 between u–v has instead two operating points, and thus each (u, k) vertices in the layered graph have two outgoing edges corresponding to the offsets k′ = ⌈k − ln Wi /∆⌉ and weights 1/Ri . This construction applies to all edges and operating points of the physical graph. QLSP finds the shortest path from (s, 0) to all (t, k) nodes, with 0 ≤ k ≤ K. In so doing, it finds the highest-rate path for all feasible end-to-end fidelities, thus the Pareto frontier.
d) Optimization: We run Dijkstra’s algorithm on the layered graph G∆ from the source (s, 0) which corresponds to the empty path with (X, Y ) = (0, 0). At the destination node t, each reachable layer k yields a candidate path pk , and the corresponding link-level operating points: the rate 1/d(t, k) is already P exact, while the exact fidelity is computed as: Y (pk ) = e∈pk ye,ae , giving Uk = U 1/d(t, k), e−Y (pk ) . The final output is the layer k ⋆ = arg maxk: d(t,k)<∞ Uk together with the path pk⋆ , and utility Uk⋆ . B. Algorithmic analysis a) Correctness: Because all layered-edge costs are nonnegative (1/Re,a ≥ 0), Dijkstra’s algorithm computes, for each layered state (v, k), the minimum cumulative inverse-rate among all paths from (s, 0) to (v, k). The layer index records the discretized cumulative Werner parameter, so the algorithm computes the best achievable rate for each Werner parameter bin and then selects the bin maximizing the e2e utility. b) Complexity: If the physical graph has n = |V| nodes and m = |E| edges, and each edge has at most A operating points, then the layered graph has |V ∆ | = O(nK) nodes and |E ∆ | = O(mKA) edges. Dijkstra’s algorithm runs in O mKA log(nK) time using a standard binary heap, and uses O(nK) memory up to storage of the layered edges. c) Approximation quality: If the cumulative Werner parameter-loss variable Y = − ln W is rounded after each edge extension using bin width ∆, then along any path of length at most L, the total discretization error in Y is at most c ≥ e−L∆ W . L∆. Hence, the e2e Werner parameter satisfies W Consider the specific case of negativity utility function (8), which, when expressed in terms of W reads U (R, W ) = R(3W − 1)/4 (we can omit the max under the assumption
2: for each e = (u, v) ∈ E do 3: for each a ∈ Ae do 4: Compute xe,a = 1/Re,a , 5: for k = 0,l. . . , K dom
ye,a = − ln We,a .
k∆+y
e,a Set k′ = . ∆ ′ 7: if k ≤ K then 8: Add layered edge (u, k) → (v, k′ ) with cost xe,a 9: if G is undirected then 10: Add layered edge (v, k) → (u, k′ ) with cost xe,a 11: Initialize distances:
6:
d(v, k) ← ∞ for all (v, k) ∈ V ∆ ,
d(s, 0) ← 0.
12: Run Dijkstra’s algorithm on G∆ from source (s, 0) 13: for each k ∈ {0, . . . , K} with d(t, k) < ∞ do 14: Recover the shortest path pk to (t, k) P 15: Compute the exact fidelity loss Y (pk ) = e∈p ye,ae . 16:
Compute Uk = U 1/d(t, k), e−Y (pk ) .
k
17: Return k ⋆ = arg maxk: d(t,k)<∞ Uk , pk⋆ , and Uk⋆ .
We2e ≥ Wmin > 1/3), the approximation ratio for the negativity term g(W ) = (3W − 1)/4 is bounded below by c) g(W 3(1 − e−L∆ ) 3L∆ ≥1− ≥1− . g(W ) 3W − 1 3Wmin − 1 The derivation uses the following inequalities: W ≤ 1 and 1− e−x ≤ x, for x ≥ 0. Thus choosing ∆ ≤ ε(3Wmin − 1)/3L is sufficient to ensure a (1 − ε)-approximation from Werner parameter discretization. For fixed Wmin bounded away from 1/3, this implies that we can set K = O (L/ε) giving a runtime of O mA Lε log nL . For simple paths, L ≤ n − 1, ε so the runtime is polynomial in the input size and 1/ε, and thus QLSP is an FPTAS for the negativity-based utility (8). The theorem below states this for general utilities. Theorem 1 (Approximation guarantee). Suppose that U is nondecreasing in rate and fidelity and that, for a given ϵ > 0, there exists some δ > 0 such that U (R, e−(Y +δ) ) ≥ (1 − ϵ)U (R, e−Y ) for all achievable (R, Y ) in the relevant domain. If QLSP uses layer width ∆ ≤ δ/L, then it returns a solution π b satisfying U (b π ) ≥ (1 − ϵ)U (π ⋆ ). The number of layers is therefore K = O(L/δ) for a fixed Wmin > 0. Thus, if 1/δ is polynomial in 1/ϵ and since L ≤ |V| − 1, QLSP is an FPTAS. Proof. The proof follows from the fact the QLSP rounds the fidelity loss of each edge upward and π ∗ uses at most L edges. Thus Yπ̂ ≤ Yπ∗ + L∆. To ensure Yπ̂ ≤ Yπ∗ + δ, it suffices to require L∆ ≤ δ. Setting ∆ = δ/L and for a fixed Wmin > 0 the number of layers K = ⌈ln(1/Wmin )/∆⌉ = O(L/δ).
To express this bound in terms of the approximation error ϵ, we must relate δ to ϵ, which depends on the form of the utility function. For the negativity utility, we derived that K = O(L/ϵ) as δ = Θ(ϵ). Similarly, for the BB84 utility, one can show that K = O (L log(1/ϵ)/ϵ). These bounds assume that Wmin is fixed and bounded away from the corresponding zeroutility threshold. Since L ≤ |V| − 1, the size of the layered graph is polynomial in K and the input size, and shortest paths can be computed in polynomial time, QLSP is an FPTAS. Fidelity-threshold utilities require additional care because an arbitrarily small decrease in fidelity can discontinuously reduce the utility to zero. This can be handled in two ways. Proposition 1 (Threshold utility). Consider Uth (R, W ) = R1[W ≥ Wreq ], where Wreq is a minimum required Werner parameter of the e2e state, and let B = − ln Wreq . Suppose an optimal solution π ⋆ satisfies Yπ⋆ ≤ B − γ for some γ > 0. If QLSP uses fidelity-loss layer width ∆ ≤ γ/L, then it returns an optimal solution for Uth . Without this margin assumption, a margin-free FPTAS can instead discretize inverse rate X = 1/R while accumulating fidelity loss Y = − ln W exactly. For any ϵ > 0, this rate-layered variant returns a feasible solution π b satisfying Yπb ≤ B and Xπb ≤ (1 + ϵ)Xπ⋆ , and hence Uth (b π ) ≥ Uth (π ⋆ )/(1 + ϵ) ≥ (1 − ϵ)Uth (π ⋆ ). C. Baselines We consider several baseline approaches. Each may suffer from a mismatch between the objective optimized during path selection and the true e2e utility, although the nature of this mismatch differs across methods. a) Sum of edge-based costs: We first define a local edge cost using one of the following choices: k • Local utility: ce = mina 1/U (Re,a , We,a ) (with k > 0), • Rate only: ce = mina 1/Re,a • Werner only: ce = − maxa ln We,a . We then minimize the sum of edge costs over the path. For local utility, the exponent k is set to the network diameter. We found this scheme to work better than simply setting k = 1 in our experiments—larger values of k emphasize fidelity over rate. After obtaining a path p, we evaluate its true e2e utility using composition rules (10). The operating point ae for each edge is determined locally based on the minimizer of the edge cost. b) Scalar surrogate sweep: An alternative approach assigns each edge the scalarized cost α − β ln We,a (13) ce (α, β) = min a∈Ae Re,a where (α, β) ≥ 0 and then runs a standard shortest-path algorithm on the graph (note: the scalar α here is not to be confused with the bright-state population parameter introduced in § III). The surrogate matches the rate and fidelity composition rules, but is based on a scalarization of the true objective. It can also be viewed as a Lagrangian relaxation of the constrained problem: min Xe2e subject to: Ye2e ≤ B, which maximizes rate subject to a fidelity constraint. The parameters
(α, β) must be swept over a range of values, and the resulting path be evaluated using the true e2e utility. The algorithm returns the path with the highest utility among all solutions found. In Theorem 2, we show that, for utilities with convex fidelity profiles, enumerating all supported solutions through this sweep recovers a globally optimal solution. We implement this sweep using an adaptive binary search. First, the coordinates X = 1/R and Y = − ln W are normalized by the median positive per-edge contributions µX and µY , respectively. Each adaptive shortest-path query therefore e + β Ye , where (X, e Ye ) = (X/µX , Y /µY ). We minimizes αX then solve the rate only and Werner only objective, corresponding to (α, β) = (1, 0) and (0, 1) respectively. Given eℓ , Yeℓ ) and zr = (X er , Yer ), two distinct solutions zℓ = (X eℓ < X er , the next shortest-path query uses the where X er − X eℓ ), for which separating weights (α, β) = (Yeℓ − Yer , X the two solutions have equal scalarized cost. If the query finds a new supported solution, the two resulting intervals are searched recursively; otherwise, that interval is exhausted. The recursion stops when no new path is found or either normalized coordinate difference is below 10−9 . Finally, every distinct feasible path identified by the sweep is evaluated using the true e2e utility, and the path with the highest utility is retained. Because the scalarized costs decompose additively over edges, for any fixed (α, β), the operating point of each selected edge can be chosen locally using (13). c) Distance path: This baseline decouples physicalroute selection from link configuration. It first selects an s– ∗ t path that P minimizes total the physical distance, pdist = arg minp e∈p Le , where Le denotes the physical length of link e. It then runs QLSP on the fixed path p∗dist using all operating points available on its links. As shown in § V, this separation is effective when the physical route can be selected largely independently of its configuration, e.g., when homogeneous links make the shortest-distance path dominant. D. Comparison Across the graphs and utility functions considered in our experiments (§ V), QLSP substantially outperforms baselines based on fixed edge costs. The surrogate sweep, however, matches QLSP when the utility function satisfies certain conditions. We formalize this result in Theorem 2, which shows that the surrogate sweep recovers a globally optimal solution when the fidelity profile is convex. Theorem 2 (Empty duality gap for convex fidelity profiles). Assume the utility has the form U (R, W ) = R · g(W ), with g nondecreasing. Define the fidelity profile as φ(y) := g e−y . Assume φ is convex and some solution achieves U > 0. Let π = (p, a) denote a physical s-t path p with ek denoting the traversed edge k on the path with a choice of operating point ak ∈ Aek . Then, on every instance the maximum of U over all routing solutions π is attained by a solution that P also minimizes the scalarized cost k Re α,a − β ln Wek ,ak k k for some weights (α, β) ≥ 0. Consequently the weightedsum sweep returns the exactly optimal utility, provided the
sweep is fine enough to hit everyP breakpoint of the piecewise linear Lagrangian ϕ(λ) = minπ k 1/Rek ,ak −λ ln Wek ,ak . There is no Lagrangian duality gap. Conversely, when φ is non-convex (e.g., g(W ) = 1[W ≥ Wreq ], where Wreq is the minimum required value, saturating ramps, sigmoids), the utility maximizer can be a Pareto-optimal solution that no scalarization ever returns, and the resulting gap can be an arbitrarily large fraction of the optimum. Proof. In the additive coordinates (X, Y ) = (1/R, − ln W ), the utility can be written as U (X, Y ) = φ(Y )/X. Convexity of φ implies that U is quasiconvex, since each sublevel set {(X, Y ) : U (X, Y ) ≤ c} = {(X, Y ) : X ≥ φ(Y )/c} is convex. Therefore, some maximizer of U over the convex hull of the achievable set lies at an extreme point. Moreover, because U is nonincreasing in both X and Y , this maximizer can be chosen on the lower-left boundary of the convex hull. Every vertex on this boundary minimizes αX + βY for some α, β ≥ 0, not both zero. Consequently, enumerating all supported solutions by weighted-sum scalarization recovers a globally optimal solution. Remark 1 (Why φ and not g). Convexity is required in the additive coordinate y = − ln W , because that is the variable in which fidelity loss accumulates linearly along a path and in which the scalarization is linear. Convexity of g in W (together with monotonicity) is sufficient but not necessary: g(W ) = √ W is concave in W yet has convex profile φ(y) = e−y/2 . Remark 2 (Instances satisfying the hypothesis). The negativity utility and the BB84 secret-key-rate utility both have convex fidelity profiles over their positive-utility regions. Consequently, for either utility, a scalarization sweep that enumerates all supported solutions recovers a globally optimal solution. Remark 3 (Frontier computation and reuse). A single run of QLSP returns the entire discretized rate–fidelity frontier d(t, k), k = 0, . . . , K, which can be re-evaluated under any nondecreasing utility and for every destination simultaneously. The scalarization sweep also admits a frontier interpretation: sweeping α, β ≥ 0 enumerates the vertices of the lower convex hull of the achievable (R, W ) set is also reusable. However, the hull omits the unsupported Pareto points (Theorem 2); it is therefore sufficient only for utilities with convex fidelity profiles, while the layered frontier remains faithful for the nonconvex ones (e.g., fidelity thresholds). The hull enumeration costs one shortest-path run per hull vertex, but the number of breakpoints of a parametric shortest path can be superpolynomial in the worst case [40], whereas the layered run is a single computation, polynomial in the input size and 1/ε. Remark 4. Our analysis assumed a discrete set of operating points A. In the Appendix, we propose a simple modification to our standard QLSP algorithm to accommodate continuous operating points; remarkably, provided a convexity condition on the per-link cost, one does not have to pay a K 2 penalty for the resulting (denser) layered graph. We show that the single-click scheme (2)-(4) satisfies the convexity condition. Furthermore, we provide in the Appendix an argument that
access to a continuous set of operating points does not eliminate the duality gap. Remark 5. Our setup and analysis assumed that s and t are allowed to choose a single network path. More generally, the nodes could multiplex between multiple network paths. In the Appendix, we show that multiplexing is not advantageous when the utility function is convex. V. N UMERICAL E VALUATION To evaluate QLSP, we employ two representative utility functions: the BB84 SKR (UBB84 ) as defined in (7), and fidelity-threshold utility (Uth ) as defined in (9), where Wreq = (4Freq − 1)/3 is the threshold value for the Werner parameter. We set Wreq = 0.94, corresponding to a Freq = 0.955, although the results are similar for other choices. These utilities illustrate both the convex-profile case (BB84) and the nonconvex case (fidelity threshold). A. Proof of Concept Fig. 2 presents a simple example of the joint decision solved by QLSP. The source s and destination t are connected by three possible corridors. The top corridor is a rate path: it is fast but comparatively noisy. The bottom corridor is a fidelity path: it is relatively low-noise but slow. The middle corridor is an adaptive path: its links offer multiple operating points (two per link; Fig. 2b), so the optimizer must determine both the physical route and the configuration of each link. The example illustrates two key effects. First, maximizing e2e utility can favor a route overlooked by optimizing rate or fidelity alone. A rate-only policy (e.g., rate only baseline from § IV-C), might select the rate path, while a fidelityonly policy (e.g., Werner only baseline from § IV-C), might select the fidelity path. Neither maximizes UBB84 . Second, even after selecting the adaptive route, link operating points cannot be chosen independently. Configuring every adaptive link for high rate produces insufficient e2e fidelity, while configuring every link for high fidelity incurs excessive delay. Instead, the optimal solution is the mixed configuration (1, 4, 1), which assigns a higher-fidelity operating point to the middle link and higher-rate operating points to the two side links. Fig. 3 uses the same network to show why baselines that use scalar edge scores can be suboptimal. For the smooth BB84 utility, the scalarized surrogate and other baselines based on local policies can recover the same adaptive configuration as QLSP. This is the favorable case: the optimum lies on a supported portion of the rate–fidelity frontier. For the threshold utility Uth with Wreq = 0.94, however, the best solution is the mixed threshold-feasible configuration (1, 4, 3): it provides just enough fidelity to satisfy the threshold while maximizing rate. The baseline instead chooses the safer fidelity path and obtains a solution worse by 22.2%. The takeaway is not that this small topology is difficult, but rather that a nonconvex utility can make the optimal configuration unsupported and therefore inaccessible to scalarization-based methods. QLSP still recovers it because it explicitly maintains a discretized e2e rate–fidelity Pareto frontier.
b) Operating points (adaptive path) 1.0
Rate path
R = 5e6, W = 0.9
r
R = 0.7e6, W = 0.995
3
a
a--b curve
sampled operating points
1 b
W link
s
0.96
t
source
4.0
s--a and b--t curves
4
0.98
Adaptive path
c) Actual end-to-end utility
continuous single-click tradeoff
Secret Key Rate / 10 5
a) Route and operating-point choice
n
available sampled point
n
QLSP-best sampled point
target
2
0.94
f
3.18 3.0 2.45
2.35 2.0
1.0
0.0
0.92
Fidelity path
3.94
3.73
0
2
4
6
R link / 10 6 pairs s -1
8
Rate path
Fid. path
Rate tune
Fid. tune
QLSP
R=2.5 W=0.81
R=0.35 W=0.99
R=1.03 W=0.884
R=0.327 W=0.965
R=0.72 W=0.927
10
1,2,1
3,4,3
1,4,1
Fig. 2. Illustration of our framework on a small network with route and operating-point choices. (a) The source s and destination t are connected by a rate-seeking path (top), a fidelity-seeking path (bottom), and an adaptive path (middle) whose links have two operating points (the two circles above each edge). (b) The operating points on the adaptive path are drawn from linear tradeoff curves, compatible with the single-click scheme. Points 1, 2 are rate-seeking points, while points 3, 4 are fidelity-seeking points; filled markers denote the best points with respect to BB84 utility. (c) End-to-end SKR utility for five candidate solutions. Rate-seeking choices lose too much fidelity, while fidelity-seeking choices sacrifice too much rate. QLSP jointly optimizes the route and operating points, selecting a mixed operating points 1, 4, 1, and achieves the largest utility. Values shown in this plot are obtained by using an approximation tolerance ϵ = 10−4 , as defined in § IV-B0c. LLEG rate and Werner parameters are calculated using the single-click scheme described in § III. a) BB84 utility
b) Fidelity-threshold utility (W > .94)
0.2
Surrogate-supported point QLSP optimum Surrogate optimum Scalarization lower hull
0.15
Shared optimum R=720k/s, W=.927
0.1
0.05 Same BB84 utility as the optimum 0.0
e2e fidelity loss Y = -ln W
e2e fidelity loss Y = -ln W
QLSP frontier 0.2
0.15
0.1 QLSP: unsupported R=450k/s, W=.946 0.05
W > 0.94
Surrogate R=350k/s, W=.990
0.0 0.5
1.0
1.5
2.0
2.5
3.0
Inverse e2e rate X = 1/R (μs/pair)
0.5
1.0
1.5
2.0
2.5
3.0
Inverse e2e rate X = 1/R (μs/pair)
Fig. 3. Limits of scalarized routing under non-convex utilities. The panels show the attainable end-to-end configurations in inverse-rate/fidelityloss space; the black dashed line is the lower convex hull reachable through surrogate scalarization. (a) For BB84 SKR, the optimum is supported, and the surrogate baseline matches QLSP. (b) Imposing the non-convex fidelity requirement W > .94 moves the optimum to an unsupported Pareto point. QLSP selects this configuration, whereas the surrogate is restricted to a suboptimal, hull-supported solution. TABLE I P ERFORMANCE BASELINES EXPRESSED AS PERCENT OF QLSP UTILITY. Topology Erdös–Rényi
Utility Local utility BB84 69.8 ± 11.2 Fidelity threshold 57.3 ± 1.1 SURFnet BB84 80.4 ± 5.5 Fidelity threshold 69.9 ± 8.1 Banded Grid BB84 0.0 ± 0.0 Fidelity threshold 0.0 ± 0.0
Distance path Werner only 86.3 ± 5.5 27.9 ± 1.5 86.8 ± 6.0 46.5 ± 2.4 92.8 ± 2.1 1.68 ± 0.06 92.7 ± 2.1 1.27 ± 0.05 100.0 ± 0.0 5.53 ± 0.08 100.0 ± 0.0 22.5 ± 3.7
B. Comparison with Baselines We evaluate the algorithms on three representative networks (Fig. 4a): a 100-node Erdös–Rényi (ER) graph [48] with p = 0.05, the 50-node SURFnet real-life topology [49], and an 8 × 12 grid with homogeneous 10 km links. For the ER topology, nodes are randomly placed on a 60 km square, and link lengths are set to the Euclidean distance between nodes. As shown in Fig. 4a, all SURFnet links use the singleclick model in (2)–(4). We set the LLEG repetition rate to Rrep = 107 s−1 , the aggregate non-fiber efficiency to c = 0.36, and the fiber loss to 0.2 dB/km, equivalently Latt ≃ 22 km. These parameters determine the rate scale but do not affect the qualitative comparisons. Each tunable single-click link is represented by A operating points sampled
uniformly over Werner parameter W ∈ [0.9, 0.9999]. We vary A to measure how performance changes as the continuous physical tradeoff is represented more finely. In the ER topology, single-click generation is fixed at W0 = 0.96 and its corresponding rate R0 . Each link instead exposes the five-point BBPSSW menu obtained from 0-4 distillation rounds using (5) and (6). The protocol-banded grid isolates hardware heterogeneity while keeping every link length equal to 10 km: links in the outer regions use the tunable single-click menu, whereas links within central columns 5–8, together with the interfaces entering and leaving this region, use the same five-point BBPSSW menu (purple in Fig. 4a). a) Discussion: The proof-of-concept in Fig. 3 established how scalarization omits solutions favored by a nonconvex utility. Fig. 4 examines whether the same phenomenon persists in representative networks. In all three settings, the QLSP frontier contains Pareto-optimal points that are not supported by its lower convex envelope. Thus, unsupported configurations arise with spatially heterogeneous links, in a real network topology, and with heterogeneous protocols. This comparison is independent of a particular utility: every displayed unsupported point is a feasible candidate that may be preferred by an appropriate monotone non-convex utility, but that no choice of surrogate weights can return. Having isolated this structural limitation of the surrogate, we compare in Table I QLSP with three heuristic baselines— Local utility, distance path, and Werner only—from § IV-C. For a baseline b, each entry reports the performance ratio ρb = 100 E(s,t) [Ub (s, t)/UQLSP (s, t)] ,
(14)
where the expectation is over all s–t pairs separated by at least eight hops. A value of 100% therefore means that the baseline always matches QLSP, while lower values measure the average fraction of QLSP utility retained. Values are obtained as an average over 50 random pairs and errors are 95% confidence intervals. We did not include rate only since it disregards fidelity entirely and its performance never matches that of QLSP.
Fig. 4. (a) Three topologies used in the evaluation. (b,c,d) The rate–fidelity frontier returned by QLSP on the three topologies. In all cases, there are frontier points (green circles) that are inside the lower convex hull and cannot thus be found by the surrogate sweep. For single click operating points, we use A = 21. a) Operating-point saturation
b) Grid-size benchmark
927.3
c) Operating-point benchmark
13.3
6.62
4.87
4.96
6.6
3.25
3.31
3.3
1.62
1.65
BB84 key rate (pairs/s)
QLSP
Solve time (s)
695.5
463.7
231.8
SURFnet
d) Tolerance benchmark
6.49
Surrogate
Distance path
10.0
Theory worst case
Banded Grid 0
0 1
5
9
13
17
Operating points per tunable link (A)
21
0 4
10
16
24
32
40
50
Grid side length (s)
1
10
25
50
75
Operating points per link (A)
100
0 10^-2
10^-3
10^-4
10^-5
10^-6
Approximation tolerance (ε)
Fig. 5. Operating-point resolution and runtime benchmarks. (a) Expected BB84 secret key rate returned by QLSP as the number of single-click operating points A increases; BBPSSW menus remain fixed. Results average 50 sampled pairs for the banded grid and 48 for SURFnet. Runtime is measured on square grids as a function of (b) grid side s, with A = 25 and ε = 10−3 ; (c) A, with s = 16 and ε = 10−3 ; and (d) ε, with s = 8 and A = 25. Runtime experiments use negativity utility. Error bars are 95% confidence intervals, with 32 timed runs per runtime point. Red curves evaluate the worst-case complexity as in § IV-B0b and are scaled as upper envelopes of the QLSP measurements. The surrogate is omitted from (d) because it does not use the layered approximation tolerance. Runtime experiments used a binary heap implementation on a MacBook Pro 14-inch (2024), 24GB.
The distance path is the strongest of these baselines, retaining approximately 86%–93% of QLSP utility on ER and SURFnet and matching QLSP on the banded grid. Selecting a physical path before configuring it can be lossless when the two decisions are separable, for example when only one path exists, when the few available paths can all be configured and compared, or when homogeneity makes a shortest path dominate independently of its configuration. The banded grid is favorable in this sense: its regular geometry makes the distance path sufficient for the sampled pairs. In general, heterogeneous link lengths or operating-point menus can make a physically longer path preferable after configuration, as the ER and SURFnet gaps demonstrate. The more local policies, as expected, are less robust: Local utility retains between 0% and approximately 80%, while Werner only retains approximately 1%–47%. Their performance also changes between BB84 and threshold utility, showing that a fixed local objective need not remain aligned with the application objective. b) Saturation and Runtime: Fig. 5a evaluates how finely the continuous single-click tradeoff must be represented. A small A can exclude useful configurations and substantially reduce the achievable key rate. SURFnet enters its high-utility regime around A = 7, whereas the banded grid stabilizes around A = 11; beyond these values, additional operating points provide only marginal gains. The resolution required for saturation depends on the topology and protocol composition. We observe similar trends across additional topologies and utility functions, omitted here for space.
Fig. 5b show that QLSP is practical at the network sizes considered here. In particular, it solves the largest 50 × 50 grid, with 2,500 nodes and 4,900 links, in approximately 11.8 seconds on a 2024 MacBook Pro. The observed growth is consistent with the polynomial worst-case analysis in § IV-B0b. Runtime grows approximately linearly with the number of operating points A (Fig. 5c), as expected from the explicit A factor in the complexity bound. Decreasing the approximation tolerance (Fig. 5d) increases the number of fidelity layers and, consequently, the explored state space; nevertheless, QLSP solves the ε = 10−6 instance on the fixed 8 × 8 grid in approximately 5.9 seconds. The theoretical curves are scaled worst-case upper envelopes rather than regression fits, so their comparison is only intended to illustrate the asymptotic trend.
The other methods reduce computational cost by searching a smaller solution space. The adaptive surrogate sweep is generally faster than QLSP, although enumerating all supported solutions may require a superpolynomial number of shortestpath queries in the worst case. The distance path restriction yields the largest reduction: on the 50 × 50 grid it requires approximately 0.28 seconds because QLSP is run only on the preselected physical path. This makes path preselection valuable when it is independent of link configuration. Outside these regimes, the runtime advantage must be weighed against the utility losses reported in Table I, since fixing the physical path can exclude the globally optimal joint route and configuration.
VI. C ONCLUSION We have presented and analyzed QLSP—an efficient approximation algorithm for jointly selecting and configuring paths in a quantum network so as to maximize a quantum utility objective between two network nodes. We have also characterized when a potentially computationally cheaper alternative scheme—based on a weighted-sum scalarization of functions of entanglement generation rate and fidelity—is optimal. Through numerical evaluation on a variety of quantum networks, we have shown that QLSP is able to reach higher quantum utility values compared to simpler baselines that track only rate, fidelity, or a local link-level utility. As future directions, QLSP could be adapted to accommodate bipartite entanglement beyond Werner states, multi-path entanglement routing (i.e., a multi-commodity formulation of the problem), and non-Pauli noise models. U SE OF AI D ISCLOSURE The authors used ChatGPT 5.6, under their supervision, to assist with documenting the code, orchestrating automated runs and collecting results, and preparing plot layouts and annotations. The authors also used ChatGPT 5.6 and Claude Opus 5 to improve the clarity of the manuscript. The authors reviewed and take responsibility for all generated content, analyses, and results. ACKNOWLEDGMENT This work was supported in part by the NSF award #2522101. It is also supported in part by the NSF grants #2346089, #2402861, and NSF- ERC Center for Quantum Networks grant EEC-1941583. R EFERENCES [1] C. H. Bennett and G. Brassard, “Quantum cryptography: Public key distribution and coin tossing,” Theoretical Computer Science, vol. 560, pp. 7–11, 2014. [2] A. K. Ekert, “Quantum cryptography based on Bell’s theorem,” Physical Review Letters, vol. 67, no. 6, pp. 661–663, aug 1991. [3] D. Gottesman, T. Jennewein, and S. Croke, “Longer-Baseline Telescopes Using Quantum Repeaters,” Physical Review Letters, vol. 109, no. 7, p. 070503, aug 2012. [4] E. T. Khabiboulline, J. Borregaard, K. De Greve, and M. D. Lukin, “Quantum-assisted telescope arrays,” Physical Review A, vol. 100, no. 2, p. 022316, 2019. [5] P. Kómár, E. M. Kessler, M. Bishof, L. Jiang, A. S. Sørensen, J. Ye, and M. D. Lukin, “A quantum network of clocks,” Nature Physics, vol. 10, no. 8, pp. 582–587, oct 2014. [6] V. Giovannetti, S. Lloyd, and L. Maccone, “Advances in quantum metrology,” Nature photonics, vol. 5, no. 4, pp. 222–229, 2011. [7] J. F. Fitzsimons and E. Kashefi, “Unconditionally verifiable blind quantum computation,” Phys. Rev. A, vol. 96, no. 1, p. 012303, 2017. [8] L. Jiang, J. M. Taylor, A. S. Sørensen, and M. D. Lukin, “Distributed quantum computation based on small quantum registers,” Physical Review A—Atomic, Molecular, and Optical Physics, vol. 76, no. 6, p. 062323, 2007. [9] G. Vardoyan and S. Wehner, “Quantum network utility maximization,” in 2023 IEEE International Conference on Quantum Computing and Engineering (QCE), vol. 1. IEEE, 2023, pp. 1238–1248. [10] C. Cabrillo, J. I. Cirac, P. Garcia-Fernandez, and P. Zoller, “Creation of entangled states of distant atoms by interference,” Physical Review A, vol. 59, no. 2, p. 1025, 1999. [11] R. F. Werner, “Quantum states with Einstein-Podolsky-Rosen correlations admitting a hidden-variable model,” Physical Review A, vol. 40, no. 8, p. 4277, 1989.
[12] M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information. Cambridge university press, 2010. [13] E. Dijkstra, “A note on two problems in connexion with graphs,” Numerische Mathematik, vol. 50, pp. 269–271, 1959. [14] G. Vidal and R. F. Werner, “Computable measure of entanglement,” Physical Review A, vol. 65, no. 3, p. 032314, 2002. [15] A. Abane, M. Cubeddu, V. S. Mai, and A. Battou, “Entanglement routing in quantum networks: A comprehensive survey,” IEEE Transactions on Quantum Engineering, 2025. [16] R. Van Meter, T. Satoh, T. D. Ladd, B. Munro, and K. Nemoto, “Path selection for quantum repeater networks,” Networking Science, vol. 3, pp. 82–95, 2013. [17] C. Di Franco and D. Ballester, “Optimal path for a quantum teleportation protocol in entangled networks,” Phys. Rev. A, vol. 85, p. 010303, Jan 2012. [18] M. Caleffi, “Optimal routing for quantum networks,” IEEE Access, vol. 5, pp. 22 299–22 312, 2017. [19] E. Schoute, L. Mancinska, T. Islam, I. Kerenidis, and S. Wehner, “Shortcuts to quantum network routing,” arXiv preprint arXiv:1610.05238, 2016. [20] M. Pant, H. Krovi, D. Towsley, L. Tassiulas, L. Jiang, P. Basu, D. Englund, and S. Guha, “Routing entanglement in the quantum internet,” npj Quantum Information, vol. 5, no. 1, p. 25, 2019. [21] S. Shi, X. Zhang, and C. Qian, “Concurrent entanglement routing for quantum networks: Model and designs,” IEEE/ACM Transactions on Networking, vol. 32, no. 3, pp. 2205–2220, 2024. [22] K. Chakraborty, D. Elkouss, B. Rijsman, and S. Wehner, “Entanglement distribution in a quantum network: A multicommodity flow-based approach,” IEEE TQE, vol. 1, pp. 1–21, 2020. [23] M. Ghaderibaneh, C. Zhan, H. Gupta, and C. Ramakrishnan, “Efficient quantum network communication using optimized entanglement swapping trees,” IEEE TQE, vol. 3, pp. 1–20, 2022. [24] G. Vardoyan, E. Van Milligen, S. Guha, S. Wehner, and D. Towsley, “On the bipartite entanglement capacity of quantum networks,” IEEE Transactions on Quantum Engineering, vol. 5, pp. 1–14, 2024. [25] J. Li, M. Wang, K. Xue, R. Li, N. Yu, Q. Sun, and J. Lu, “Fidelityguaranteed entanglement routing in quantum networks,” IEEE Transactions on Communications, vol. 70, no. 10, pp. 6748–6763, 2022. [26] Y. Zhao, G. Zhao, and C. Qiao, “E2E fidelity aware routing and purification for throughput maximization in quantum networks,” in IEEE INFOCOM 2022, 2022, pp. 480–489. [27] Z. Xiao, J. Li, K. Xue, N. Yu, R. Li, Q. Sun, and J. Lu, “Purification scheduling control for throughput maximization in quantum networks,” Communications Physics, vol. 7, no. 1, p. 307, 2024. [28] M. Victora, S. Tserkis, S. Krastanov, A. S. de la Cerda, S. Willis, and P. Narang, “Entanglement purification on quantum networks,” Physical Review Research, vol. 5, no. 3, p. 033171, 2023. [29] M. Skrzypczyk and S. Wehner, “An architecture for meeting quality-ofservice requirements in multi-user quantum networks,” 2021. [Online]. Available: https://arxiv.org/abs/2111.13124 [30] F. P. Kelly, K. Aman, Maulloo, and D. K. H. Tan, “Rate control for communication networks: shadow prices, proportional fairness and stability.” Journal of the Operational Research Society, vol. 49, no. 3, pp. 237–252, 1998. [31] S. H. Low and D. E. Lapsley, “Optimization flow control. I. Basic algorithm and convergence.” IEEE/ACM Transactions on Networking,, vol. 7, no. 6, pp. 861–874, 1999. [32] D. P. Palomar and M. Chiang, “Alternative distributed algorithms for network utility maximization: Framework and applications.” IEEE Trans. on Automatic Control, vol. 52, no. 12, pp. 2254–2269, 2007. [33] N. K. Panigrahy, L. Bacciottini, C. Hollot, E. A. Van Milligen, M. G. de Andrade, N. S. Rao, G. Vardoyan, and D. Towsley, “A framework for distributed resource allocation in quantum networks,” arXiv preprint arXiv:2510.09371, 2025. [34] S. Kar and A. Mukhopadhyay, “On Utility-optimal Entanglement Routing in Quantum Networks,” in 2026 International Conference on Quantum Communications, Networking, and Computing (QCNC). IEEE, 2026, pp. 457–464. [35] Q. Zhang, N. Di Cicco, M. Ibrahimi, R. C. Almeida Jr., A. Gatto, R. Boutaba, and M. Tornatore, “Link configuration for fidelityconstrained entanglement routing in quantum networks,” in IEEE INFOCOM 2025 – IEEE Conference on Computer Communications, 2025.
[36] Z. Wang and J. Crowcroft, “Quality-of-service routing for supporting multimedia applications,” IEEE Journal on Selected Areas in Communications, vol. 14, no. 7, pp. 1228–1234, 1996. [37] R. Hassin, “Approximation schemes for the restricted shortest path problem,” Mathematics of Operations Research, vol. 17, no. 1, pp. 36– 42, 1992. [38] D. H. Lorenz and D. Raz, “A simple efficient approximation scheme for the restricted shortest path problem,” Operations Research Letters, vol. 28, no. 5, pp. 213–219, 2001. [39] A. Jüttner, B. Szviatovszki, I. Mécs, and Z. Rajkó, “Lagrange relaxation based method for the QoS routing problem,” in Proc. IEEE INFOCOM, 2001, pp. 859–868. [40] P. J. Carstensen, “The complexity of some problems in parametric linear and combinatorial programming,” Ph.D. dissertation, University of Michigan, 1983. [41] B. C. Coutinho, R. Monteiro, L. Bugalho, and F. A. Monteiro, “Entanglement routing based on fidelity curves,” arXiv preprint arXiv:2303.12864, 2024. [42] W. J. Munro, K. Azuma, K. Tamaki, and K. Nemoto, “Inside quantum repeaters,” IEEE Journal of Selected Topics in Quantum Electronics, vol. 21, no. 3, pp. 78–90, 2015. [43] C. H. Bennett, G. Brassard, S. Popescu, B. Schumacher, J. A. Smolin, and W. K. Wootters, “Purification of noisy entanglement and faithful teleportation via noisy channels,” Physical review letters, vol. 76, no. 5, p. 722, 1996. [44] D. Deutsch, A. Ekert, R. Jozsa, C. Macchiavello, S. Popescu, and A. Sanpera, “Quantum privacy amplification and the security of quantum cryptography over noisy channels,” Physical review letters, vol. 77, no. 13, p. 2818, 1996. [45] C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K. Wootters, “Mixed-state entanglement and quantum error correction,” Phys. Rev. A, vol. 54, pp. 3824–3851, Nov 1996. [46] P. G. Kwiat, E. Waks, A. G. White, I. Appelbaum, and P. H. Eberhard, “Ultrabright source of polarization-entangled photons,” Physical Review A, vol. 60, no. 2, p. R773, 1999. [47] W. Dür, H.-J. Briegel, J. Cirac, and P. Zoller, “Quantum repeaters based on entanglement purification,” Phys. Rev. A, vol. 59, no. 1, p. 169, 1999. [48] P. Erdős and A. Rényi, “On random graphs I,” Publicationes Mathematicae Debrecen, vol. 6, pp. 290–297, 1959. [49] S. Knight, H. X. Nguyen, N. Falkner, R. Bowden, and M. Roughan, “The internet topology zoo,” IEEE Journal on Selected Areas in Communications, vol. 29, no. 9, pp. 1765–1775, 2011. [50] A. Aggarwal, M. Klawe, S. Moran, P. Shor, and R. Wilber, “Geometric applications of a matrix searching algorithm,” in Proceedings of the second annual symposium on Computational geometry, 1986, pp. 285– 292.
A PPENDIX
T IME MULTIPLEXING Generally, the network may time-multiplex among multiple routes or configurations. Consider two routing solutions with end-to-end rate–fidelity pairs (R1 , F1 ) and (R2 , F2 ). Suppose that the first solution is used for a fraction f ∈ [0, 1] of the time. The aggregate generation rate is R(f ) = f R1 + (1 − f )R2 .
Because the two solutions generate f R1 and (1 − f )R2 pairs per unit time, respectively, the fidelity of the average distributed state is F (f ) =
f R1 F1 + (1 − f )R2 F2 . f R1 + (1 − f )R2
(16)
Consider a utility of the form U (R, F ) = R g(F ), where g is convex. Let λ(f ) = f R1 /(f R1 + (1 − f )R2 ). By Jensen’s inequality, U (R(f ), F (f )) = R(f )g(λ(f )F1 + (1 − λ(f ))F2 ) ≤ R(f ) [λ(f )g(F1 ) + (1 − λ(f ))g(F2 )] = f R1 g(F1 ) + (1 − f )R2 g(F2 ) = f U (R1 , F1 ) + (1 − f )U (R2 , F2 ) ≤ max{U (R1 , F1 ), U (R2 , F2 )}.
(17)
Thus, time multiplexing cannot outperform the better of the two solutions when g is convex. This includes the negativity and BB84 utilities. Time multiplexing can, however, improve a nonconvex utility. Consider Uth (R, F ) = R1[F ≥ Freq ]. Suppose solution 1 has higher rate but falls below the threshold, while solution 2 is feasible: R1 > R2 ,
F1 < Freq ≤ F2 .
The multiplexed solution satisfies the fidelity requirement if and only if f R1 (F1 − Freq ) + (1 − f )R2 (F2 − Freq ) ≥ 0.
In this appendix, we show that: 1) time multiplexing between multiple network paths is not advantageous over using exclusively one path, when the utility is convex; 2) continuous operating points do not eliminate the duality gap, so although scalarization allows us to incorporate continuous operating points easily, it is not guaranteed to recover the optimum; and 3) under certain conditions, incorporating continuous operating points into the layered graph can be done without paying the K 2 penalty. The basic idea is to switch to a hop-indexed dynamic program (DP), similar to Bellman–Ford, and exploit a special property of the transition structure to compute all of the min-plus values efficiently. When the transition matrix has the so-called Monge structure, one can show that the indices of the optimal transitions vary monotonically across layers. This allows all optimal values to be computed theoretically in O(K) time, or in O(K log K) time using a simpler divideand-conquer scheme.
(15)
(18)
Since R(f ) increases with f , the optimal mixture uses the largest feasible fraction of the high-rate solution. The threshold is therefore active, yielding f⋆ =
R2 (F2 − Freq ) . R1 (Freq − F1 ) + R2 (F2 − Freq )
(19)
The resulting rate is ⋆ Rmux =
R1 R2 (F2 − F1 ) . R1 (Freq − F1 ) + R2 (F2 − Freq )
(20)
Whenever F2 > Freq , we have f ⋆ > 0, and hence ⋆ Rmux = f ⋆ R1 + (1 − f ⋆ )R2 > R2 .
Thus, mixing a high-rate, sub-threshold solution with a lowerrate, high-fidelity solution can strictly outperform the feasible solution used alone.
C ONTINUOUS O PERATING P OINTS Continuous operating points do not eliminate the global scalarization gap. Even when the configuration problem on every fixed path is convex and has zero duality gap, the full routing problem includes a discrete choice among physical paths. The global achievable set is a union of path-specific feasible sets and need not be convex. Moreover, each path may require a different Lagrange multiplier. Thus, a configuration can be optimal for the constrained problem yet unsupported by every global weighted-sum scalarization. Fig. 6 illustrates this reasoning. The network in this example contains three internally disjoint two-link paths P1 , P2 , and P3 between s and t. Solutions are shown from coarse to relatively fine operating point tuning (as the plots go from left to right). In the (X, Y ) = (1/R, − ln W ) plane, the configurations of each physical path form a cluster of points. With continuous link tuning, these clusters extend into path-specific attainable regions. Such regions can overlap; they happen to be disjoint here. Increasing the resolution of operating points fills these regions more densely, but does not make their union convex. All links follow the single-click model (2)–(4), with Rrep = 108 s−1 . The effective transmissions on the first links of P1 , P2 , and P3 are, respectively, (η1 , η2 , η3 ) = 10−3 (1, 0.8, 10/3); the second link of Pi has transmission 4ηi . Both links on each path are independently tunable over W ∈ [0.900, 0.910], [0.925, 0.940], and [0.988, 0.992], respectively. These pathdependent values are chosen for illustrative purposes so that the attainable regions of these paths do not overlap. We sample A = 3, 17, 257 uniformly spaced operating points per link, and use Wmin = 0.8 and K = 4463 fidelity-loss bins. Utilities are evaluated using the recovered e2e Werner parameter; exhaustive enumeration provides the sampled frontier and checks the QLSP and adaptive-surrogate solutions. For Uth with Wreq = 0.88, QLSP selects an unsupported configuration (teal stars in the figure) on P2 at every sampled resolution, whereas the surrogate selects P3 . At A = 257, their rates are approximately 6.05 × 103 and 4.80 × 103 pairs/s, respectively. Finer sampling thus resolves the intermediate island but does not bring its sampled optimum onto the lower convex envelope. E FFICIENT LAYERED GRAPHS APPROACH FOR CONTINUOUS OPERATING POINTS
Continuous operating points Ae can be incorporated without discretizing the operating-point curve when the per-link inverse-rate cost has suitable convex structure in y = − ln W . In this case, each edge relaxation is a convex min-plus convolution whose transition matrix is Monge. A hop-indexed dynamic program can compute all layer transitions through an edge in O(K log K) time using a simple divide-andconquer algorithm, or in O(K) time using the “SMAWK” algorithm [50]. The resulting complexities are O(L|E|K log K) and O(L|E|K), respectively, rather than the O(L|E|K 2 ) cost of explicitly examining all pairs of layers.
Let Dh (v, k) denote the minimum cumulative inverse rate X = 1/R of any path from s to v that uses at most h physical edges and terminates in fidelity-loss layer k. For each edge e ∈ E, define 1 , (21) ce (j) = min a∈Ae : Re,a ⌈− ln We,a /∆⌉=j
with ce (j) = ∞ if no operating point induces layer increment j. Here Ae may be a continuous feasible set of operating points. The dynamic program is initialized as D0 (s, 0) = 0,
D0 (v, k) = ∞
for (v, k) ̸= (s, 0). (22)
For h = 1, . . . , L, and for k = 0, . . . , K, define the relaxation through an edge e = (u, v) as Th,e (k) = min {Dh−1 (u, k − j) + ce (j)} . 0≤j≤k
(23)
The equation above is a min-plus convolution. The dynamicprogramming recurrence is then Dh (v, k) = min Dh−1 (v, k), min Th,e (k) . (24) e=(u,v)∈E
After L iterations, each reachable destination layer k yields a candidate path with rate 1 . (25) Rk = DL (t, k) Its layer represents fidelity loss approximately k∆, or Wk ≈ e−k∆ . As in the standard QLSP algorithm, the exact fidelity of the recovered path can instead be computed from its selected operating points and used to evaluate the final utility. When ce (j) is convex in j, the edge relaxation (23) is a Monge min-plus convolution, allowing all destination-layer values to be computed jointly rather than examining all O(K 2 ) source–destination layer pairs. a) Convexity of the per-link cost: The Monge acceleration applies whenever the minimum inverse-rate cost ce (j) is convex in the fidelity-loss increment j. This condition holds, in particular, for the continuous single-click operating-point model in Eqs. (2)–(4). Recall that 4 F (α) = 1 − α, R(α) = 2Rrep ηe α, W (α) = 1 − α, 3 where ηe = η(Le ) is the transmissivity of edge e. Eliminating α gives the affine rate–Werner tradeoff 3 Re (W ) = Rrep ηe (1 − W ). (26) 2 Writing y = − ln W , so that W = e−y , the inverse-rate cost becomes 1 2 xe (y) = = . (27) Re (y) 3Rrep ηe (1 − e−y ) This function is decreasing and strictly convex (in the relevant region of y, which is nonnegative for W ∈ [0, 1]): 2e−y < 0, 3Rrep ηe (1 − e−y )2 2e−y (1 + e−y ) x′′e (y) = > 0. 3Rrep ηe (1 − e−y )3 x′e (y) = −
(28) (29)
e2e fidelity loss Y = -ln W
0.20
a) Coarse tuning
b) Fine tuning
c) Near-continuous tuning
A = 3 operating points per link
A = 17 operating points per link
A = 257 operating points per link
1
P1
s
t
2
P2 0.15
W ≥ 0.88
0.10
P3
0.05 0.00
0.20
0.20
0.15
0.15
Unsupported even with continuous tuning
3
0.10 0.15 0.20 0.25 0.30 Inverse e2e rate X = 1/R (ms/pair)
Attainable configurations QLSP optimum
W ≥ 0.88
0.10 0.05 0.00
W ≥ 0.88
0.10 0.05
0.10 0.15 0.20 0.25 0.30 Inverse e2e rate X = 1/R (ms/pair)
Sampled Pareto frontier Best surrogate solution
0.00
0.10 0.15 0.20 0.25 0.30 Inverse e2e rate X = 1/R (ms/pair)
Scalarization-supported Lower convex envelope (bound)
Fig. 6. Refining operating points within fixed physical paths. The same three paths are sampled with a) A = 3, b) A = 17, and c) A = 257 operating points per link. Pale points show attainable configurations; hollow circles mark sampled Pareto points and orange circles mark scalarization-supported points. The horizontal line marks Wreq = 0.88, with feasible configurations below it. Stars identify QLSP solutions; larger orange circles identify the best feasible surrogate solutions.
Since xe (y) is decreasing, the minimum inverse-rate cost within fidelity-loss bin j is attained at the largest feasible value of y in that bin. For bins whose upper endpoint j∆ lies within the feasible operating range, ce (j) = xe (j∆).
(30)
Sampling a convex function on a uniform grid preserves discrete convexity, and hence 2ce (j) ≤ ce (j − 1) + ce (j + 1).
(31)
b) Monge structure: For a fixed edge e = (u, v) and hop iteration h, define the transition matrix over feasible source– destination layer pairs as Mh,e (k, i) = Dh−1 (u, i) + ce (k − i).
i
(33)
Mh,e (k, i) + Mh,e (k + 1, i + 1) (34)
so the transition matrix has Monge structure over its feasible entries. Consequently, if i⋆ (k) denotes the smallest sourcelayer index attaining the minimum for destination layer k, then i⋆ (k) ≤ i⋆ (k + 1).
im = i⋆ (km ). Monotonicity implies that all destination layers to the left of km have minimizers in [iℓ , im ], while all layers to the right have minimizers in [im , ir ]. We therefore recurse on [kℓ , km − 1] × [iℓ , im ] and
which is obtained by re-parameterizing j = k − i in (23). Then, discrete convexity in (31) gives
≤ Mh,e (k, i + 1) + Mh,e (k + 1, i),
by scanning the feasible source indices in [iℓ , ir ] and obtaining
(32)
The edge relaxation is its row minimum: Th,e (k) = min Mh,e (k, i),
c) An O(K log K) divide-and-conquer algorithm: The monotonicity property (35) yields a simple divide-and-conquer algorithm for computing all row minima. Suppose we wish to compute Th,e (k) for destination layers k ∈ [kℓ , kr ] and know that their minimizing source layers lie in i ∈ [iℓ , ir ]. We first evaluate the middle layer kℓ + kr km = 2
(35)
Thus, the optimal source-layer index moves monotonically forward as the destination layer increases.
[km + 1, kr ] × [im , ir ]. At each recursion depth, the total number of candidate source indices examined is O(K), and there are O(log K) depths. Hence all transitions through one physical edge can be computed in O(K log K) time using this simple implementation, giving an overall hop-indexed complexity of O(L|E|K log K) .
(36)
The SMAWK linear-time Monge matrix-search algorithm [50] further reduces each edge relaxation to O(K) and the overall complexity to O(L|E|K).