1
Inter-Satellite Link Optimization for Low-Latency Global Networking
arXiv:2604.15528v1 [cs.NI] 16 Apr 2026
Arman Mollakhani, Jerayu Tiamraj, Shu-Jie Cao, and Dongning Guo
Abstract—Large-scale low-Earth-orbit satellite constellations offer a promising platform for global low-latency networking, aided by faster propagation in free space than in fiber and copper. In such systems, end-to-end latency is largely determined by the inter-satellite link (ISL) topology. In particular, the network diameter—the maximum shortest path between any pair of satellites—serves as a key performance metric for time-sensitive applications. Designing diameter-optimal topologies is challenging due to degree constraints, line-of-sight limitations, and orbital dynamics. This paper proposes a two-stage optimization framework for ISL topology design. First, a continuous relaxation of the link selection problem is formulated as a convex program that maximizes the algebraic connectivity of the Laplacian, serving as a tractable surrogate for diameter minimization. Second, the resulting fractional solution is mapped to a feasible discrete topology using integer linear programming. An iterative local-search heuristic is also developed as a baseline. Extensive simulations on Walker–Delta constellations show that the proposed method consistently achieves smaller network diameters and improved robustness compared to conventional heuristics, while allowing trade-offs between latency and link persistence. The approach offers a principled framework for designing high-performance satellite mesh networks. For a constellation of 1,500 satellites, each equipped with four ISLs of up to 2,500 km, the network diameter can be reduced to as low as 12, yielding end-to-end delays under 90 ms between any two points on Earth.
I. I NTRODUCTION Low-Earth-orbit (LEO) constellations such as Starlink, AST SpaceMobile, OneWeb, and Amazon Leo aim to provide nearubiquitous global coverage [1]. Optical inter-satellite links (ISLs) enable efficient in-network data transport within and across constellations. In free space, electromagnetic waves propagate at approximately 300 km per millisecond (ms), traversing the Earth’s pole-to-pole geodesic distance in as little as 67 ms—roughly 47% faster than in optical fiber [2]. Consequently, a well-designed LEO ISL network can outperform the fastest terrestrial fiber routes in end-to-end latency, making such systems well-suited for latency-critical applications including remote control, telemedicine, high-frequency trading, blockchain synchronization, and emergency communications [3], [4]. Constellations operated by a single service provider can further exploit globally optimized routing in stead of distributed protocols such as the Border Gateway Protocol (BGP): forwarding tables can be precomputed for short epochs, dissemThe authors are with the Department of Electrical and Computer Engineering, Northwestern University, Evanston, IL 60208. Emails: {arman.mollakhani, shujie.cao, dguo}@northwestern.edu, [email protected] This work was supported in part by NSF grants Nos. 2132700 (SpectrumX) and 2434044.
inated to all nodes, and updated as the topology evolves [2]. This deterministic approach eliminates in-orbit route discovery and enables techniques such as cut-through switching [5] or multiprotocol label switching forwarding [6], reducing perhop processing overhead to the microsecond scale [7]. In this work, we use the network diameter—the maximum shortest-path cost between any two satellites—as a structural proxy for latency. While the framework accommodates arbitrary positive link costs, we adopt hop count as the path metric in our simulations. This choice is justified by the geometric structure of LEO constellations and the bounded length of ISLs. Consequently, paths with fewer hops tend to be more direct, reducing propagation distance. In practice, each satellite is limited to a small number of steerable laser terminals, and any candidate link must satisfy distance and line-of-sight constraints not just at a single instant but ideally over an entire orbital period. The design problem is therefore to select a sparse, physically feasible set of ISLs that achieves a small diameter while maintaining orbital stability. Most existing approaches rely on fixed geometric patterns or on iterative local-search heuristics. These methods often fail to capture the global structural properties required for near-optimal connectivity. Moreover, local-search methods are prone to poor local minima. This paper addresses these limitations through a spectral graph-theoretic framework. The main contributions of this paper are as follows: (i) We formulate a constrained network diameter minimization problem that accounts for hardware limits, geometric constraints, and long-term link stability. (ii) We propose a two-stage spectral optimization framework that combines convex relaxation and integer linear programming to design degree-constrained ISL topologies. (iii) We develop an efficient first-order projected gradient ascent method that scales to large constellations. (iv) We introduce an iterative local-search heuristic that serves both as a practical standalone method and as a warm start for the spectral framework. The remainder of this paper is organized as follows. Sec. II reviews related work. Sec. III describes the constellation model. Sec. IV formulates the diameter-minimization problem. Sec. V introduces the spectral framework. Sec. VI presents the heuristic algorithm. Sec. VII provides simulation results, and Sec. VIII concludes the paper. II. R ELATED W ORK Prior work on ISL topology design has explored timedivision methods for periodically activating feasible links for
2
long-term connectivity [8], link scheduling for navigation constellations [9], and multi-objective assignment algorithms that account for payload and visibility constraints. Large-scale ISL scheduling has also been modeled as a discrete network multicommodity flow problem, with data-driven search heuristics proposed to minimize transmission delay [10]. While these approaches emphasize temporal reachability, they do not explicitly optimize worst-case latency or network diameter. Multicast and broadcast strategies in LEO constellations have also been explored, including core-based multicast trees and disruption-tolerant mechanisms for multilayered satellite networks [11], [12]. Surveys of routing in LEO networks highlight trade-offs between centralized and distributed protocols and the need to adapt to predictable topological changes [13]. Centralized, software-defined networking approaches precompute routing tables for anticipated topologies, reducing onboard computation and routing convergence time [14]. While improving resilience and group communication, these methods often rely on idealized routing assumptions and ignore geometric constraints inherent to LEO dynamics. Stochastic geometry has been used to model spatial connectivity and evaluate average-case latency [15], and congestionaware routing protocols improve throughput under dynamic traffic [16]. Analytic models estimate ISL hop counts between ground users [17], and algorithms have been proposed to compute shortest-distance paths directly from satellite phase [18]. In [19], we have proposed an iterative local-search heuristics to refine inter-plane link assignments in Starlink-inspired constellations, but the method is prone to local optima as constellation size grows [20]. Spectral graph theory provides a framework to relate local link selection to global network properties. The algebraic connectivity of the Laplacian captures network expansion and diameter [21], [22], and maximizing this spectral gap via convex optimization can eliminate bottlenecks [23]. Continuous relaxations of discrete link selection problems, supported by semidefinite programming (SDP) and convex analysis [24]– [26], enable tractable optimization of low-diameter networks. We leverage these spectral techniques to design ISL topologies under practical constraints. III. S YSTEM M ODEL In this work, we model the continuous evolution of an LEO satellite network using an epoch-based framework, where the topology is assumed to remain static within each epoch, and is only updated at discrete time intervals. In other words, each epoch represents a fixed configuration of feasible ISLs, during which routing and link assignment decisions can be made deterministically. This formulation allows for tractable analysis and static optimization of network properties within each epoch.
θ, giving orbital radius r = RE + h, where RE = 6,371 km is the mean radius of Earth. Orbital planes are evenly spaced in longitude, and satellites within each plane are uniformly distributed. To capture practical position variations, a random phase offset ϕi ∈ [0, ϕmax ] is applied to each plane i. This preserves the overall structure while modeling deviations due to orbit insertion errors and gravitational perturbations that cause drift from the nominal geometry over time [27]. B. Link Feasibility Constraints With the origin of a three-dimensional Cartesian space placed at Earth’s center, we denote the position of satellite u by xu ∈ R3 . In the time-varying case, the position at time t is denoted xu (t). The maximum feasible distance for an ISL is denoted by dmax , determined by physical-layer constraints including laser power, beam divergence, and receiver sensitivity. A candidate link between two satellites u and v is considered feasible at time t if the following conditions are satisfied: 1) Distance constraint: The Euclidean distance must not exceed the maximum link distance: ∥xu (t) − xv (t)∥ ≤ dmax .
(1)
2) Line-of-sight constraint: The link must not be obstructed by the Earth. Mathematically, this requires ∥xu (t) × xv (t)∥ > RE , ∥xu (t) − xv (t)∥
(2)
where the left-hand side of (2) computes the perpendicular distance from Earth’s center to the line segment connecting u and v (since the vector cross product xu (t) × xv (t) gives a vector whose length is equal to the area of the parallelogram formed by those two vectors). C. Snapshot vs. Viability-Constrained Models We require each selected ISL to remain feasible over an epoch, i.e., geometrically feasible, free of Earth obstruction, and within distance bounds for the entire epoch. Accordingly, we adopt an epoch-based model in which the routing topology remains static within each epoch, allowing route precomputation and reducing in-network overhead. We consider two different link feasibility models: Definition 1 (Snapshot Model). A link between two satellites u and v is feasible if it satisfies constraints (1) and (2) at a single instant in time t0 . The set of all such feasible links is denoted E.
A. Constellation Geometry
Definition 2 (Viability-Constrained Model). A link between two satellites u and v is feasible only if it satisfies constraints (1) and (2) for the entire duration of an orbital period T , i.e., for all t ∈ [t0 , t0 + T ]. The resulting feasible set, denoted E ′ , satisfies E ′ ⊆ E.
We consider a Walker–Delta constellation with Np orbital planes and Ns satellites per plane, for a total of N = Np ×Ns satellites. Each satellite orbits at altitude h km with inclination
This distinction is practically important because reconfiguring optical ISLs is slow. Retargeting requires physically steering laser terminals, re-acquiring beam alignment, with
3
setup times on the order of several seconds [28]. These delays limit frequent reconfiguration, favoring topologies that remain valid over longer durations to maintain stable and low-latency performance. IV. P ROBLEM F ORMULATION Let V be fixed and denote the set of all N satellites. Let F denote the set of all feasible (undirected) ISLs and E ⊂ F denote the set of active ISLs. We represent the satellite network as an undirected graph G = (V, E). Each satellite is assumed to be equipped with D high-speed, bidirectional optical terminals for inter-satellite communication. Consequently, the total degree of each satellite is bounded: degE (v) ≤ D,
∀v ∈ V.
(3)
The general problem is to select a subset of edges E ⊂ F that optimizes some performance metric subject to degree constraints (3). Here, we let the metric be the worst-case network diameter. To define this precisely, each edge (u, v) ∈ F is assigned a fixed cost pu,v ∈ (0, +∞). For all (u, v) ∈ / F, it is convenient to assume pu,v = +∞. We define a path between satellites u and v in graph (V, E) as a sequence of vertices P = (v1 , v2 , . . . , vm ), where v1 = u, vm = v, and (vi , vi+1 ) ∈ E for i = 1, . . . , m − 1. . The total cost of the path P is: C(P ) =
m−1 X
pvi ,vi+1 .
(4)
Snapshot Model: With F = E, (7) minimizes the diameter subject to geometric and line-of-sight constraints (1) and (2) evaluated at a single instant t0 . ′ • Viability-Constrained Model: With F = E , (7) minimizes the diameter subject to the stricter requirement that selected links satisfy the geometric and line-of-sight constraints for an entire orbital period T . Finding a globally optimal solution to (7) is computationally prohibitive for large-scale constellations in general. This follows because even the simpler problem of finding a minimum-diameter degree-constrained spanning subgraph is known to be NP-hard [20]. To see this, note that the decision version of our problem—“does there exist a subgraph satisfying all constraints with diameter at most D?”—subsumes degree-constrained subgraph feasibility as a special case. We pursue two complementary approaches to solve (7): (i) a principled spectral optimization framework that provides continuous relaxation and discrete rounding (Sec. V), and (ii) an effective heuristic algorithm based on iterative local search (Sec. VI). Throughout this paper, we set pu,v = 1 for all edges, so that path costs reduce to hop counts. •
V. A T OPOLOGY D ESIGN F RAMEWORK To overcome the computational intractability of the discrete diameter minimization problem (7), we propose a two-stage graph-theoretic framework that leverages continuous relaxation and spectral surrogate optimization.
i=1
Let Pu,v (E) denote the set of all paths between u and v. The shortest-path cost between u and v, which depends on E, is defined as: dE (u, v) =
min
P ∈Pu,v (E)
C(P ).
(5)
If no path exists between u and v, then Pu,v (E) is empty and dE (u, v) = +∞. Definition 3 (Network Diameter). Given the set of edges E, the network diameter is defined as: diam(E) = max dE (u, v). u,v∈V
(6)
Given a universal vertex degree constraint D, the discrete optimization problem seeks to select a subset of edges E ⊂ F that minimizes the network diameter: minimize
diam(E)
subject to
degE (v) ≤ D,
E⊂F
(7a) ∀v ∈ V.
(7b)
The edge cost can be modeled in different ways: • Weighted graph model: Each edge in E is assigned a positive weight representing the latency, which may include propagation, transmission, and queuing delays. • Hop-count model: Each edge is assigned one unit of cost, making dE (u, v) the shortest-path hop distance between u and v. The candidate set F takes different forms depending on the system’s operational requirements (as outlined in Sec. III):
A. Spectral Optimization and Continuous Relaxation We first introduce a continuous relaxation by associating a strength variable xij ∈ [0, 1] with each potential edge (i, j) ∈ F . Since the graph is undirected, xij = xji . The degree constraint is relaxed to a continuous budget: X xij ≤ D, ∀i ∈ V. (8) j:(i,j)∈F
Next, rather than directly minimizing the worst-case path cost, we employ a spectral surrogate objective to maximize the graph’s overall connectivity. We construct the N ×N weighted Laplacian matrix L(X), where the symmetric weight matrix X contains the link strengths xij : (P if i = j k̸=i xik Lij = (9) −xij if i ̸= j. The connection between minimizing graph diameter and maximizing the second-smallest eigenvalue of the Laplacian, λ2 (the algebraic connectivity or spectral gap), is rigorously grounded in spectral graph theory. A graph with a large λ2 exhibits strong expansion properties, indicating that no sparse cut exists to separate the graph into two large components. This is formalized by Cheeger’s inequality [21], which bounds the Cheeger constant h(G): p λ2 ≤ h(G) ≤ 2λ2 . (10) 2 For families of graphs with bounded degrees, such expansion properties are known to imply logarithmic upper bounds on
4
the graph diameter [21], [22]. Therefore, maximizing λ2 serves as a principled, mathematically tractable proxy for eliminating bottlenecks and minimizing the overall network diameter. We therefore formulate the diameter minimization problem as the following spectral optimization problem: maximize {xij }
subject to
λ2 (L(X)) X xij ≤ D,
(11a) ∀i ∈ V
diagonal matrix of edge weights. The weighted Laplacian matrix L satisfies: L = BΛB T .
Proof. Let A = BΛB T , and let Aij be the (i, j)-th entry of this matrix. Since Λ is diagonal, Aij =
(11b)
xij = xji ,
∀(i, j) ∈ F ∀(i, j) ∈ F.
M X
Λkk Bik Bjk .
(14)
k=1
j:(i,j)∈F
0 ≤ xij ≤ 1,
(13)
(11c) (11d)
A high strength xij → 1 increases the link’s contribution to overall graph connectivity. The spectral problem (11) can be solved efficiently (e.g., by being cast as an SDP). Proposition 1 (Convexity). The spectral optimization problem (11) is a convex optimization problem. Proof. Let the eigenvalues of L be λ1 ≤ λ2 ≤ · · · ≤ λN . Since λ1 = 0 for any connected graph, maximizing λ2 is equivalent to minimizing (−λ1 − λ2 ), which is the sum of the two largest eigenvalues of the matrix −L. The sum of the k largest eigenvalues of a symmetric matrix is a wellknown convex function of that matrix’s entries [24], [29], [30]. Therefore, the objective is a convex function of xij ’s. The constraints are linear and hence form a convex set. B. Efficient Sparse Formulation While problem (11) is convex, a direct implementation using an N × N matrix variable involves O(N 2 ) decision variables. For large constellations, this would be computationally prohibitive, rendering standard SDP solvers prohibitively memory-intensive. However, the graph is naturally sparse; the set of feasible links F is much smaller than the set of all pairs. To exploit this sparsity, we reformulate the problem using a vector of variables. Let M = |F |. We index the potential edges k = 1, . . . , M . We define the decision variable as a vector x ∈ RM , where xk corresponds to the strength of the k-th potential edge connecting nodes i and j. Furthermore, let N (i) denote the set of indices of edges incident to vertex i in F . To construct the Laplacian efficiently, we utilize the incidence matrix B ∈ RN ×M . For each edge k connecting nodes u and v, we assign an arbitrary direction (e.g., u → v) such that: if i = u 1 Bik = −1 if i = v (12) 0 otherwise. We define the diagonal weight matrix Λ ∈ RM ×M as Λ = diag(x). We can now state the following relationship between the vector of edge weights and the Laplacian matrix. Lemma 1 (Laplacian Decomposition). Let B be the unweighted incidence matrix of a graph, and let Λ be the
We analyze two cases: 1) Off-diagonal terms (i ̸= j): The term Bik Bjk is non-zero only if edge k connects vertices i and j. If it does, one entry is 1 and the other is −1, yielding a product of −1. Thus, Aij = −xij if edge (i, j) exists, and 0 otherwise. This matches the definition of Lij for i ̸= j. PM 2) Diagonal terms (i = j): Here, Aii = k=1 Λkk (Bik )2 . Since Bik ∈ {0, 1, −1}, (Bik )2 = 1 if edge k is connected to node i, and 0 otherwise. Thus, Aii = P Λ , which is the weighted degree of node i. kk k∈N (i) This matches the definition of Lii . Therefore, A = L. Using this lemma, we express L as a linear function of the vector x: L(x) = B diag(x)B T .
(15)
This formulation reduces the number of optimization variables from O(N 2 ) to O(M ), making the problem tractable for large constellations. C. First-Order Optimization While the spectral problem (11) is convex, its SDP formulation requires enforcing positive semidefiniteness of an N × N matrix at each solver iteration. Standard primal-dual interior-point methods for SDPs scale with complexity O(N 3 ) to O(N 4 ) per iteration [26]. To address this, we propose a firstorder projected gradient ascent (PGA) method [31]. PGA is well-suited for this problem because: i) it operates directly on the sparse vector x, avoiding the O(N 2 ) memory footprint of SDPs; and ii) the projection onto the box constraints 0 ≤ xij ≤ 1 is computationally trivial. We reformulate the constrained optimization problem as a sequence of maximizations using a quadratic penalty method. 1) Penalty Formulation: To relax the hard degree constraints (8), we introduce a quadratic penalty function. Recall that N (i) denotes the set of neighbors of node i in the initial graph (V, F ). Denote the current total strength at node i as X si (x) = xij , (16) j∈N (i)
and the penalty as Φ(x) =
X
2
(max (0, gi (x))) ,
(17)
i∈V
where gi (x) = si (x) − D.
(18)
5
We then seek to maximize a parameterized Lagrangian-like objective. The penalized objective, for scalar penalty parameter ρ > 0, is: maximize x∈[0,1]M
J (x; ρ) = f (x) − ρ Φ(x),
(19)
where f (x) = λ2 (L(x)).
(20)
3) Penalty Gradient: We now derive the gradient of the penalty term defined in (17). For a specific edge e = (u, v), which contributes to the degree sums of both node u and node v, we have ∇xuv Φ(x) ∂ X 2 (max(0, gk (x))) = ∂xuv
(26)
k∈V
As established in penalty method theory [31], [32], as ρ → ∞, the solution sequence x∗ (ρ) converges to the solution of the original constrained problem. Intuitively, as ρ grows, any nonzero violation of the constraints incurs an arbitrarily high cost, forcing the optimizer into the feasible region to maximize the net objective. 2) Spectral Gradients and Eigenvalue Multiplicity: The function defined by (20) is concave but can be non-smooth when λ2 has multiplicity greater than one. In highly symmetric topologies such as Walker constellations, eigenvalues frequently coalesce, such that λ2 = λ3 = · · · = λk . At such points, the standard gradient is undefined because the eigenvectors are not unique; they form an invariant subspace, and the function exhibits a non-differentiable “kink” [25], [33]. To proceed, we utilize the concept of the generalized gradient (or Clarke subdifferential) for spectral functions [34]. Let ϵ be a numerical tolerance. We define the cluster of active eigenvalues K as: K = {k ∈ {2, . . . , N } : |λk − λ2 | < ϵ}.
(21)
We first consider the derivative for a simple eigenvalue. Let v(k) be the normalized eigenvector associated with eigenvalue λk . From standard matrix perturbation theory [35], the differential of a simple eigenvalue λk with respect to the matrix L is given by dλk = (v(k) )T (dL)v(k) . Recall that L = B diag(x)B T , where x is the vector of edge weights. The partial derivative with respect to a specific edge weight xij connecting nodes i and j is: ∂λk ∂L = (v(k) )T v(k) (22) ∂xij ∂xij = (v(k) )T (ei − ej )(ei − ej )T v(k) (23) 2 (k) (k) (24) = vi − v j where ei is the standard basis vector. When the multiplicity is greater than 1, the subdifferential ∂f (x) is defined as the convex hull of the gradients induced by all unit vectors in the eigenspace of λ2 [25]. To define a robust ascent direction, we approximate the subgradient by averaging the gradients over the active cluster K. This technique, a standard heuristic in eigenvalue optimization [33], stabilizes the trajectory by accounting for the sensitivity of the entire subspace: 2 1 X (k) (k) vi − v j . (25) ∇xe f (x) ≈ |K| k∈K
Averaging produces a direction in the convex hull that is invariant to arbitrary rotations of the eigenbasis within the cluster, preventing the optimizer from oscillating wildly between different eigenvectors.
∂gu ∂gv + 2 max(0, gv (x)) ∂xuv ∂xuv = 2 max(0, su (x) − D) + 2 max(0, sv (x) − D)
= 2 max(0, gu (x))
(27) (28)
∂gk is equal to 1 if k = u where we have used the fact that ∂x uv or k = v and 0 otherwise. The net gradient direction gt at iteration t combines the spectral pull and the penalty push:
gt = ∇f (xt ) − ρt ∇Φ(xt ).
(29)
4) Momentum-Based Update and Annealing Schedule: To escape local optima and traverse the flat plateaus characteristic of spectral functions, we employ a momentum-based update rule, specifically Polyak’s heavy-ball method [36]. Momentum accelerates convergence by accumulating velocity in directions of persistent descent and dampening oscillations in directions of high curvature [37]. We simultaneously employ a dynamic annealing schedule. Let Tmax be the total number of iterations. We define a variable penalty parameter ρt that ramps quadratically to strictly enforce feasibility in the final iterations, and a learning rate ηt that decays to ensure the variance of the gradient approximation vanishes: 2 t ρt = ρmin + (ρmax − ρmin ) (30) Tmax η0 ηt = . (31) 1 + αt Let µ ∈ [0, 1) denote the momentum coefficient, and let mt denote the momentum vector, initialized to the zero vector. Let P[0,1] denote the Euclidean projection onto the unit hypercube, defined element-wise as min(1, max(0, x)). The complete update rule is: mt+1 = µmt + ηt gt
(32)
xt+1 = P[0,1] (xt + mt+1 ) .
(33)
D. The Two-Stage Solution Methodology Our proposed solution combines the continuous spectral relaxation with a final discrete selection step. a) Stage 1: Spectral Relaxation.: We solve the spectral optimization problem (11) using the PGA method described above. Let {x∗ij } be the optimal strengths. b) Stage 2: Discrete Rounding via integer linear programming (ILP): We use these optimal x∗ij ’s as the weights in our discrete optimization. We introduce binary decision variables yij , where yij = 1 signifies that edge (i, j) is selected. The problem becomes to select a set of edges that
6
maximizes the total spectral strength under the original degree budget: X maximize x∗ij · yij (34a) {yij }
subject to
(i,j)∈F
X
yij ≤ D,
∀i ∈ V
(34b)
∀(i, j) ∈ F
(34c)
j:(i,j)∈F
yij ∈ {0, 1}, yij = yji ,
∀(i, j) ∈ F.
(34d)
This two-stage approach is more tractable. It replaces a single, highly complex discrete problem with two sequential, betterunderstood problems: a convex optimization (Stage 1) and a combinatorial optimization (Stage 2) for which efficient solvers and heuristics exist.
E. Specialization: Fixed Intra-Plane Links We classify links as intra-plane, connecting satellites within the same orbital plane, and inter-plane, connecting satellites across different planes. While all links in F are treated equally in our optimization framework and no links are pre-assigned, it also accommodates configurations where specific links are fixed by the hardware architecture. In certain operational scenarios, intra-plane links may be mandated by hardware design or mission requirements. For example, in constellations where each satellite carries dedicated laser terminals for its two orbital-plane neighbors (as in Starlink [38], [39]), these links are always active and not subject to optimization. The framework readily accommodates this as a special case. Specifically, let Eintra and Einter denote the intra-plane and inter-plane subsets of F , respectively. Intra-plane links are fixed with strength xij = 1, and the optimization is performed only over Einter , subject to a separate per-satellite inter-plane budget Dinter = D − 2 (accounting for the two fixed intraplane links). Mathematically, this is achieved by appending an additional constraint to the spectral optimization problem (11): xij = 1,
∀(i, j) ∈ Eintra .
(35)
Stage 1 then solves the spectral relaxation over the inter-plane variables only, with the fixed intra-plane contributions included as constants in the Laplacian. Stage 2 rounds the resulting weights via the ILP with budget Dinter . This specialization reduces the number of decision variables from |F | to |Einter | and can be employed whenever the intra-plane topology is predetermined by the constellation’s hardware architecture.
VI. A N I TERATIVE L OCAL -S EARCH A LGORITHM We now describe an iterative local-search heuristic approach [19] that serves as a baseline for evaluating the spectral framework. The heuristic operates over the same candidate edge set F defined in Sec. IV, with identical feasibility constraints and degree budget D.
A. Initial Topology Construction The algorithm constructs a baseline network graph (V, E0 ) where E0 ⊂ F using the following procedure. We begin with the set of vertices V with no links. The topology is constructed via a greedy, degree-balanced strategy. The following four-step sequence is executed iteratively over multiple passes to ensure maximum utilization of the degree budget D: 1) All satellites with degree less than D are identified and sorted by residual capacity (available link ports), so that satellites with fewer existing links are prioritized for new connections. 2) For a selected satellite u, all feasible neighbors v ∈ F such that degE (v) < D are identified as candidates. 3) To promote long-range connections that reduce network diameter, a stratified selection strategy is employed. Candidate neighbors are sorted by Euclidean distance and split into “far” and “near” halves. Each half is independently shuffled to promote diversity. The algorithm then attempts to add links starting with the far set, followed by the near set, until u reaches its degree budget or no candidates remain. 4) A link (u, v) is established only if both u and v have not yet reached their respective link limits and the link does not already exist. B. Iterative Topology Optimization After initialization, an iterative local-search algorithm refines the ISL connections to reduce the network diameter diam(E) as defined in (6). The optimization alternates between two phases over a fixed number of n iterations. 1) Repair phase (every K iterations): Every K-th iteration, a repair phase is triggered to address connectivity gaps. The algorithm identifies all “deficient” satellites with degree less than D. For each such satellite, it attempts to add new links by connecting to valid candidates that have available capacity. To ensure fairness, the candidate list is shuffled before attempting connections. A safeguard mechanism ensures that no satellite exceeds its degree constraint as a result of this process. This phase corrects structural weaknesses and prevents sparsely connected nodes from becoming communication bottlenecks. 2) Random replacement phase (default iterations): In all non-repair iterations, the algorithm performs a random replacement step to explore the solution space. A random subset of m satellites is selected, and for each, one of its existing links (if any) is randomly removed. The algorithm then attempts to form a new link with a different, randomly selected valid candidate, respecting all constraints. This phase introduces topological perturbations that allow the search to escape local minima and discover more efficient configurations. 3) Evaluation and acceptance: After each iteration, the modified topology E ′ is evaluated and compared against the best-found topology E ∗ using an ordered multi-objective criterion. Let H̄ denotes the average maximum shortest-path cost over all satellites, and σ denotes the stable ISL fraction. The new topology E ′ replaces E ∗ if: (a) It has a smaller diameter: diam(E ′ ) < diam(E ∗ ); or
7
Algorithm 1 Iterative Link Optimization [19] 1: Input: Satellite set V , feasible link set F , degree budget D, total iterations n, reinforcement interval K, perturbation size m 2: Output: Optimized topology E ∗
for the spectral framework: the discrete topology produced by Algorithm 1 is converted into an initial strength vector x0 that initializes the PGA of Sec. V-C, ensuring that the spectral method improves upon or matches the heuristic baseline.
3: E ∗ ← initial topology according to Section VI-A 4: Evaluate E ∗ : diameter H ∗ , avg max hop H̄ ∗ , stability σ ∗ 5: for i = 1 to n do 6: E ′ ← an identical copy of E ∗ 7: if i mod K = 0 then ▷ Repair Phase 8: U ← {v ∈ V : degE ′ (v) < D}
VII. S IMULATION AND R ESULTS
for each satellite u ∈ U do Add feasible links to valid candidates Enforce degree constraints at both endpoints end for else ▷ Random Replacement Phase for each u in a random subset of m satellites do Replace one random link (u, v) (if any) by feasible link (u, w) end for end if Evaluate E ′ : diameter H ′ , avg max hop H̄ ′ , stability
9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19:
σ′ ′
∗
′
∗
′
∗
′
if H < H or (H = H and H̄ < H̄ ) or (H = H ∗ and H̄ ′ = H̄ ∗ and σ ′ > σ ∗ ) then 21: E ∗ ← E ′ , H ∗ ← H ′ , H̄ ∗ ← H̄ ′ , σ ∗ ← σ ′ 22: end if 23: end for 24: return E ∗ 20:
(b) It has the same diameter but a smaller average shortestpath cost: diam(E ′ ) = diam(E ∗ ) and H̄ ′ < H̄ ∗ ; or (c) Both primary metrics are equal, but the fraction of links that remain stable over an orbital period is higher: diam(E ′ ) = diam(E ∗ ), H̄ ′ = H̄ ∗ , and σ ′ > σ ∗ . This multi-objective acceptance criterion ensures greedy improvement on the primary objective (diameter reduction) while using secondary metrics for tie-breaking to enhance overall network performance and stability. C. Operational Variants The baseline algorithm is guided by two variants with distinct link feasibility rules. In the snapshot variant, the set of valid candidates for each satellite is larger, as links need only be feasible at a single instant. This greater topological flexibility allows link stability (σ) to be used as a secondary tie-breaking metric. In contrast, the viability-constrained variant only permits links that remain geometrically valid over a full orbital period. In this stricter model, link stability is guaranteed by construction and therefore does not serve as a separate optimization criterion. The complete baseline procedure is summarized in Algorithm 1. This adaptive local-search method navigates the solution space to find topologies balancing low network diameter with long-term link stability. While effective as a standalone method, the baseline also serves as a warm-start
We evaluate the performance of the proposed two-stage spectral framework against the baseline heuristic. A. Simulation Setup a) Constellation Parameters.: Simulations are conducted on Walker–Delta constellations where each orbit has an independent phase offset uniformly distributed in [0, ϕmax ]. The parameters for a 500-satellite constellation and a 1,584satellite constellation are summarized in Table I. The baseline heuristic runs for n = 300 iterations per trial, with a repair phase every K = 15 iterations and a perturbation size of m = 20 satellites. TABLE I C ONSTELLATION AND S YSTEM PARAMETERS Parameter
Small-Scale
Large-Scale
Total satellites (N ) Orbital planes (Np ) Satellites per plane (Ns ) Orbital altitude (h) Orbital inclination (θ) Total degree budget (D) Maximum ISL distance (dmax ) Maximum orbital phase offset (ϕmax )
500 25 20 550 km 53◦ 4 3,500 km π/4
1,584 72 22 550 km 53◦ 4 Varied π/4
b) Solver Configuration.: The spectral framework’s PGA is configured with a maximum of Tmax = 20,000 iterations. The dynamic penalty schedule (30) initializes with ρmin = 0.001 to allow temporary constraint violations during the early exploration phase, and gradually increases (anneals) to ρmax = 30.0 to enforce feasibility more strictly in later steps. The learning rate (31) uses η0 = 2.0 with decay α = 0.002, and the momentum coefficient is µ = 0.8. Additionally, the tolerance ϵ used to identify the active eigenvalue cluster K in (21) is adapted dynamically at each step t as ϵt = 0.05|λ2 (xt )| + 10−4 . c) Warm-Start Strategy.: The baseline heuristic is first executed to obtain a feasible initial topology. This discrete solution is converted into an initial strength vector x0 for the PGA. d) Statistical Methodology.: All experiments are repeated over 50 independent random initializations. We report minimum, maximum, and mean values to characterize statistical variance. We report three metrics for each scenario/algorithm: • the maximum hop count (diameter), • the average maximum hop count over all satellites (as a source), and • the total number of edges in the resulting subgraph. Each topology is evaluated using all-source shortest-path analysis via breadth-first search (BFS) [40].
8
B. Small-Scale Validation (N = 500)
Spectral (min–max) Heuristic (mean)
15
Heuristic (min–max)
14 Hops
Table II compares the methods in the snapshot scenario. The proposed spectral framework consistently achieves a worstcase diameter of 9 hops across all 50 runs, whereas the heuristic baseline described in Sec. VI averages 11.02 hops. Furthermore, the spectral method consistently producing a 4regular graph with 1,000 edges, whereas the heuristic method leaves some vertices with fewer than four edges.
Spectral (mean)
16
13 12 11
TABLE II S NAPSHOT R ESULTS (N = 500 ) Spectral Framework
Heuristic Baseline
Diameter (Hops) Min Max Mean
9.0 9.0 9.0
11.0 12.0 11.02
Avg. Max Hops Min Max Mean
8.11 8.20 8.16
9.81 10.42 10.07
Total Edges Min Max Mean
1000.0 1000.0 1000.0
990.0 999.0 995.38
2,000
2,200
2,400
2,600
2,800
3,000
Max ISL Distance (km) Fig. 1. Avg. maximum hops vs. ISL range (snapshot scenario, N = 1,584). Solid lines show the mean over 50 runs; shaded bands indicate the min–max range.
20
Table III presents results under viability constraints, where E ′ is restricted to links stable over a full orbital period. The spectral framework achieves a consistent diameter of 11, compared to 13.08 for the heuristic.
Spectral (mean)
19
Hops
Metric
10
Spectral (min–max)
18
Heuristic (mean)
17
Heuristic (min–max)
16 15 14 13
TABLE III V IABILITY-C ONSTRAINED R ESULTS (N = 500)
12 11
Metric
Spectral Framework
Heuristic Baseline
Diameter (Hops) Min Max Mean
11.0 11.0 11.0
13.0 14.0 13.08
Avg. Max Hops Min Max Mean
11.0 11.0 11.0
11.82 12.37 12.08
Total Edges Min Max Mean
1000.0 1000.0 1000.0
987.0 999.0 993.42
2,000
2,200
2,400
2,600
2,800
3,000
Max ISL Distance (km) Fig. 2. Avg. maximum hops vs. ISL range (viability-constrained scenario, N = 1,584). Solid lines show the mean over 50 runs; shaded bands indicate the min–max range.
18
Spectral (mean)
17
Spectral (min–max) Heuristic (mean)
C. Large-Scale Scalability Analysis (N = 1,584) To assess scalability, we simulate a dense constellation with 72 orbital planes and 22 satellites per plane (N = 1,584), with each plane assigned a random phase shift. The maximum ISL range dmax is varied from 2,000 km to 3,000 km. The results are summarized in Figures 1–4, which plot the minimum, maximum, and mean values across 50 runs. To translate hop counts into latency, consider the snapshot case with dmax = 2,500 km as an example: a mean diameter of 12 hops yields end-to-end delays under 90 ms, since most hops are below the maximum range.
Hops
16
Heuristic (min–max)
15 14 13 12 11 10 2,000
2,200
2,400
2,600
2,800
3,000
Max ISL Distance (km) Fig. 3. Diameter vs. ISL range (snapshot scenario, N = 1,584). Solid lines show the mean over 50 runs; shaded bands indicate the min–max range.
Hops
9
22 21 20 19 18 17 16 15 14 13 12
Spectral (mean) Spectral (min–max) Heuristic (mean) Heuristic (min–max)
leaves satellites under-connected. This is a structural advantage of the continuous relaxation, which distributes strength smoothly across edges before the rounding step ensures a regular graph. VIII. C ONCLUSION
2,000
2,200
2,400
2,600
2,800
3,000
Max ISL Distance (km) Fig. 4. Diameter vs. ISL range (viability-constrained scenario, N = 1,584). Solid lines show the mean over 50 runs; shaded bands indicate the min–max range.
We have presented a computationally efficient framework for designing sparse, low-diameter ISL topologies in largescale LEO constellations. Simulations on Starlink-like Walker– Delta constellations show that the spectral framework consistently reduces network diameter by 2–3 hops relative to the heuristic approach. Future directions include incorporating distance-weighted path costs to directly target propagation delay, adaptive topologies responding to dynamic traffic, and multi-shell architectures with relay satellites for further diameter reduction. ACKNOWLEDGMENTS
D. Discussion The simulation results yield several key observations. a) Spectral advantage: Across all configurations, the spectral framework consistently achieves lower network diameters than the heuristic baseline. At N = 500 in the snapshot scenario, the spectral method achieves a diameter of 9, while the heuristic averages 11.02—a reduction of approximately 18%. At N = 1,584 with dmax = 2,500 km, the spectral method achieves a mean diameter of 12.08 compared to 14.98 for the heuristic in the snapshot case, and 15.00 vs. 17.02 in the viability-constrained case. This improvement is attributable to the spectral method’s ability to optimize for global graph properties (expansion, bottleneck elimination) rather than relying on local perturbations. b) Stability–latency trade-off: The results clearly illustrate the fundamental tension between latency performance and operational stability. The snapshot model yields lowerdiameter networks by exploiting transient, short-term links. However, this comes at the cost of high link churn, which would necessitate frequent topology updates and complex dynamic routing protocols, potentially negating the latency benefits. In contrast, the viability-constrained model produces completely stable topologies with higher diameters. Such static or semi-static topologies are well-suited for centrally managed mega-constellations where routing tables can be pre-computed and link reconfiguration is minimized. c) Effect of ISL range: Increasing dmax improves both diameter and λ2 under both scenarios. Each 100 km increase in dmax yields roughly a 0.3–0.5 hop reduction in mean diameter at the large scale. This suggests that improvements in laser terminal technology (enabling longer-range links) would directly translate to latency improvements. The gains are more pronounced in the viability-constrained case, where longerrange links are more likely to remain stable across orbital periods. d) Degree budget utilization: The spectral framework consistently achieves full degree budget utilization (all satellites maintain exactly D links), while the heuristic occasionally
The authors thank Dr. Aravindan Vijayaraghavan for suggesting the spectral graph-theoretic approach. R EFERENCES [1] R. Berry, P. Bustamante, D. Guo, T. Hazlett, M. Honig, I. Murtazashvili, S. Palo, and M. H. Weiss, “Spectrum rights in outer space: Interference management for low earth orbit (LEO) broadband constellations,” Journal of Information Policy, vol. 14, 2024. [2] M. Handley, “Delay is not an option: Low latency routing in space,” in Proceedings of the 17th ACM Workshop on Hot Topics in Networks, 2018, pp. 85–91. [3] I. N. Bozkurt, A. Aguirre, B. Chandrasekaran, P. B. Godfrey, G. Laughlin, B. Maggs, and A. Singla, “Why is the internet so slow?!” in International Conference on Passive and Active Network Measurement. Springer, 2017, pp. 173–187. [4] A. Mollakhani and D. Guo, “Fault-tolerant spectrum usage consensus for low-earth-orbit satellite constellations,” in 2025 IEEE International Conference on Decentralized Applications and Infrastructures (DAPPS). IEEE, 2025, pp. 21–26. [5] P. Kermani and L. Kleinrock, “Virtual cut-through: A new computer communication switching technique,” Computer Networks (1976), vol. 3, no. 4, pp. 267–286, 1979. [6] E. Rosen, A. Viswanathan, and R. Callon, “Multiprotocol label switching architecture,” Internet Engineering Task Force (IETF), RFC 3031, Jan 2001. [7] S. Heine, C. Hofmann, and A. Knopp, “End-to-end latency evaluation of LEO satellite IoT systems,” in IET Conference Proceedings CP873, vol. 2023, no. 48. IET, 2023, pp. 216–221. [8] X. Chu and Y. Chen, “Time division inter-satellite link topology generation problem: Modeling and solution,” International Journal of Satellite Communications and Networking, vol. 36, no. 2, pp. 194–206, 2018. [9] L. Sun, J. Yang, W. Huang, L. Xu, S. Cao, and H. Shao, “Inter-satellite time synchronization and ranging link assignment for autonomous navigation satellite constellations,” Advances in Space Research, vol. 69, no. 6, pp. 2421–2432, 2022. [10] J. Liu, L. Xing, L. Wang, Y. Du, J. Yan, and Y. Chen, “A datadriven parallel adaptive large neighborhood search algorithm for a largescale inter-satellite link scheduling problem,” Swarm and Evolutionary Computation, vol. 74, p. 101124, 2022. [11] L. Cheng, J. Zhang, and K. Liu, “Core-based shared tree multicast routing algorithms for LEO satellite IP networks,” Chinese Journal of Aeronautics, vol. 20, no. 4, pp. 353–361, 2007. [12] G. Zheng and Y. Guo, “A routing strategy with link disruption tolerance for multilayered satellite networks,” International Journal of Communications, Network and System Sciences, vol. 3, no. 11, pp. 835–842, 2010. [13] X. Qi, J. Ma, D. Wu, L. Liu, and S. Hu, “A survey of routing techniques for satellite networks,” Journal of Communications and Information Networks, vol. 1, no. 4, pp. 66–85, 2016.
10
[14] M. Corici, H. Buhr, H. Zope, and M. Zaboub, “An SDN-based solution for mega-constellation routing,” in 2024 IEEE 35th International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC). IEEE, 2024, pp. 1–6. [15] R. Wang, M. A. Kishk, and M.-S. Alouini, “Stochastic geometrybased low latency routing in massive LEO satellite networks,” IEEE Transactions on Aerospace and Electronic Systems, vol. 58, no. 5, pp. 3881–3894, 2022. [16] S. Dai, L. Rui, S. Chen, and X. Qiu, “A distributed congestion control routing protocol based on traffic classification in LEO satellite networks,” in 2021 IFIP/IEEE International Symposium on Integrated Network Management (IM). IEEE, 2021, pp. 523–529. [17] Q. Chen, G. Giambene, L. Yang, C. Fan, and X. Chen, “Analysis of inter-satellite link paths for LEO mega-constellation networks,” IEEE Transactions on Vehicular Technology, vol. 70, no. 3, pp. 2743–2755, 2021. [18] Q. Chen, L. Yang, Y. Zhao, Y. Wang, H. Zhou, and X. Chen, “Shortest path in LEO satellite constellation networks: An explicit analytic approach,” IEEE Journal on Selected Areas in Communications, vol. 42, no. 5, pp. 1175–1187, 2024. [19] A. Mollakhani, J. Tiamraj, S.-J. Cao, and D. Guo, “Inter-Satellite Link Configuration for Fast Delivery in Low-Earth-Orbit Constellations,” in Proceedings of the 2026 IEEE Aerospace Conference (AeroConf). Big Sky, MT, USA: IEEE, Mar. 2026. [20] M. R. Garey and D. S. Johnson, Computers and Intractability: A Guide to the Theory of NP-Completeness. San Francisco, CA, USA: W. H. Freeman, 1979. [21] F. R. K. Chung, Spectral Graph Theory. American Mathematical Society, 1997, vol. 92. [22] ——, “Diameters and eigenvalues,” Journal of the American Mathematical Society, vol. 2, no. 2, pp. 187–196, 1989. [23] A. Ghosh and S. Boyd, “Growing well-connected graphs,” in Proceedings of the 45th IEEE Conference on Decision and Control. IEEE, 2006, pp. 6605–6611. [24] K. Fan, “On a theorem of weyl concerning eigenvalues of linear transformations i,” Proceedings of the National Academy of Sciences, vol. 35, no. 11, pp. 652–655, 1949. [25] A. S. Lewis, “Convex analysis on the hermitian matrices,” SIAM Journal on Optimization, vol. 6, no. 1, pp. 164–177, 1996. [26] L. Vandenberghe and S. Boyd, “Semidefinite programming,” SIAM Review, vol. 38, no. 1, pp. 49–95, 1996. [27] M. Hu, F. Li, W. Xue, C. Liu, W. Guo, and Y. Ruan, “Station maintenance for low-orbit large-scale constellations based on absolute and relative control strategies,” Applied Sciences, vol. 15, no. 9, p. 4640, 2025. [28] D. Bhattacharjee, A. U. Chaudhry, H. Yanikomeroglu, P. Hu, and G. Lamontagne, “Laser inter-satellite link setup delay: Quantification, impact, and tolerable value,” in 2023 IEEE Wireless Communications and Networking Conference (WCNC). IEEE, 2023, pp. 1–6. [29] A. S. Lewis, “The convex analysis of unitarily invariant matrix functions,” Journal of Convex Analysis, vol. 2, no. 1/2, pp. 173–183, 1996. [30] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge University Press, 2004. [31] D. P. Bertsekas, Nonlinear Programming. Belmont, MA: Athena Scientific, 1999. [32] J. Nocedal and S. J. Wright, Numerical Optimization. Springer, 1999. [33] M. L. Overton and R. S. Womersley, “On the sum of the largest eigenvalues of a symmetric matrix,” SIAM Journal on Matrix Analysis and Applications, vol. 13, no. 1, pp. 41–45, 1992. [34] F. H. Clarke, Optimization and Nonsmooth Analysis. SIAM, 1990. [35] G. W. Stewart and J.-g. Sun, Matrix Perturbation Theory. Academic Press, 1990. [36] B. T. Polyak, “Some methods of speeding up the convergence of iteration methods,” USSR Computational Mathematics and Mathematical Physics, vol. 4, no. 5, pp. 1–17, 1964. [37] N. Qian, “On the momentum term in gradient descent learning algorithms,” Neural Networks, vol. 12, no. 1, pp. 145–151, 1999. [38] Federal Communications Commission, “Request for Modification of the Authorization for the SpaceX NGSO Satellite System,” Report and Order, FCC 21-48A1, 2021. [Online]. Available: https: //docs.fcc.gov/public/attachments/FCC-21-48A1.pdf [39] C.-J. Wang, “Structural properties of a low Earth orbit satellite constellation—the Walker Delta network,” in Proceedings of MILCOM’93—IEEE Military Communications Conference, vol. 3. IEEE, 1993, pp. 968–972. [40] D. C. Kozen, “Depth-first and breadth-first search,” in The Design and Analysis of Algorithms. Springer, 1992, pp. 19–24.