Paper deep dive
Neural operator learning for collision-aware trajectory planning of spacecraft swarms
Sidhdharth D. Sikka, Suyi Gao, Zehui Lu, Rongjie Lai, Shaoshuai Mou
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 8/4/2026, 3:22:45 AM
Summary
This paper introduces a permutation-equivariant neural operator for collision-aware trajectory planning of spacecraft swarms in congested orbits. The method maps distributions of spacecraft, targets, and debris to collision-free trajectories in a single forward pass, trained via self-supervised physics objectives and adversarial threats. It generalizes zero-shot to swarms of 1,000 agents amid over 11,000 catalogued objects, matching optimal-control accuracy while significantly reducing computational cost compared to classical optimization methods.
Entities (8)
Relation Signals (6)
Permutation-equivariant neural operator → solves → collision-aware trajectory planning
confidence 95% · Here we introduce a permutation-equivariant neural operator that maps distributions of spacecraft, targets and debris to collision-aware trajectories
Permutation-equivariant neural operator → outperforms → debris-blind baseline
confidence 93% · evading worst-case threats that a debris-blind baseline cannot
Permutation-equivariant neural operator → generalizesto → swarms of 1,000 agents
confidence 92% · Trained on ten spacecraft, it generalizes zero-shot to swarms of 1,000
Gauss-Newton finish → refines → neural operator output
confidence 90% · The operator is trained... paired with a batched Gauss-Newton finish that enforces exact orbital dynamics
Two-Line Element (TLE) catalog → provides → initial states
confidence 88% · every scenario draws its NN spacecraft initial states from the real Two-Line Element (TLE) catalog
SpaceX Starlink → demonstrates → high collision-avoidance maneuver rate
confidence 85% · SpaceX’s Starlink constellation executed more than 207,000 automated collision-avoidance maneuvers
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Autonomous spacecraft swarms must plan fuel-efficient, collision-free maneuvers in increasingly congested orbits, yet classical trajectory optimization scales poorly as pairwise safety constraints multiply with swarm size, and learning-based planners rarely transfer across swarm sizes or debris densities. Here we introduce a permutation-equivariant neural operator that maps distributions of spacecraft, targets and debris to collision-aware trajectories for an entire swarm in a single forward pass, paired with a batched Gauss-Newton finish that enforces exact orbital dynamics. The operator is trained without optimal-trajectory labels, combining self-supervised physics objectives with adversarial threats generated against its own rollouts. Trained on ten spacecraft, it generalizes zero-shot to swarms of 1,000 amid more than 11,000 catalogued objects, matching a per-agent optimal-control solver's accuracy, evading worst-case threats that a debris-blind baseline cannot, and reducing proximity within the swarm several-fold. Physics-grounded operator learning thus offers a fast, scalable alternative to optimal control for crowded orbits.
Tags
Links
- Source: https://arxiv.org/abs/2608.00320v1
- Canonical: https://arxiv.org/abs/2608.00320v1
Trouble viewing inline? Open PDF directly →
Full Text
75,032 characters extracted from source content.
Expand or collapse full text
[1,2] D. [1] of Aeronautics and Astronautics, University, Lafayette, , 2] Research Group, Francisco, , 3] of Mathematics, University, Lafayette, , 4] Researcher, , , Neural operator learning for collision-aware trajectory planning of spacecraft swarms sikkas@purdue.edu gao757@purdue.edu zehuilu789@gmail.com lairj@purdue.edu mous@purdue.edu * [ [ [ Abstract Autonomous spacecraft swarms must plan fuel-efficient, collision-free maneuvers in increasingly congested orbits, yet classical trajectory optimization scales poorly as pairwise safety constraints multiply with swarm size, and learning-based planners rarely transfer across swarm sizes or debris densities. Here we introduce a permutation-equivariant neural operator that maps distributions of spacecraft, targets and debris to collision-aware trajectories for an entire swarm in a single forward pass, paired with a batched Gauss–Newton finish that enforces exact orbital dynamics. The operator is trained without optimal-trajectory labels, combining self-supervised physics objectives with adversarial threats generated against its own rollouts. Trained on ten spacecraft, it generalizes zero-shot to swarms of 1,000 amid more than 11,000 catalogued objects, matching a per-agent optimal-control solver’s accuracy, evading worst-case threats that a debris-blind baseline cannot, and reducing proximity within the swarm several-fold. Physics-grounded operator learning thus offers a fast, scalable alternative to optimal control for crowded orbits. keywords: Neural operators, Operator learning, Distributional control, Multi-agent systems, Collision-aware planning, Spacecraft swarms, Orbital autonomy Low Earth Orbit (LEO) is becoming increasingly crowded as satellite constellations expand and debris accumulates. Roughly 16,000 active satellites and an estimated 140 million debris fragments now occupy near-Earth space [1]. Between December 2025 and May 2026 alone, SpaceX’s Starlink constellation executed more than 207,000 automated collision-avoidance maneuvers, over three times its rate a year earlier, averaging more than 40 maneuvers per satellite per year, nearly one per week [2]. As orbital congestion increases, trajectory planning becomes a persistent, large-scale coordination problem rather than an occasional corrective maneuver. Frequent avoidance maneuvers reduce mission lifetime, consume fuel, and create coordination burdens that grow with constellation size. Classical approaches to multi-agent orbital maneuvering rely primarily on centralized optimization: mixed-integer linear programs offer safety guarantees [3], linearized and distributed convex optimization improves tractability [4], and model predictive and receding-horizon schemes handle constraints in closed loop [5, 6], but all must re-solve programs whose collision constraints multiply with the number of agents and debris objects. Reactive schemes such as velocity obstacles [7], artificial potential fields [8], and deployed rule-based screening avoid this cost but grow conservative and prone to mutual conflict in dense traffic [2], while heuristic [9], reinforcement-learning [10], and decision-theoretic [11] planners, including learned swarm-navigation policies [12], are flexible but seldom transfer across swarm sizes and debris densities. Figure 1: Neural operator framework for collision-aware swarm planning. The spacecraft swarm and debris field are represented as distributions over orbital elements. Separate attention pathways encode target alignment and debris avoidance, while a physics-informed baseline (green dashed) plus learned residual (blue solid) generates dynamically grounded, collision-aware trajectories for the full swarm. As a final step, the resulting debris-avoiding trajectories are refined by a batched Gauss–Newton finish that closes each one onto exact two-body dynamics. A key difficulty is structural. Collision constraints couple each spacecraft to other spacecraft and debris objects through pairwise interactions, causing the number of safety relationships to grow rapidly with swarm size and debris density. Rather than treating a swarm as a collection of individual spacecraft, one can instead model it as a probability distribution evolving under controlled dynamics. Mean Field Games (MFGs) formalize this perspective by characterizing the limit of infinitely many interacting agents through coupled Hamilton–Jacobi–Bellman and Fokker–Planck equations [13]. While this representation scales more naturally with population size and has shown promise in large-scale robotic systems [14], solving the associated partial differential equations remains computationally demanding in high-dimensional physical domains, even with dedicated machine-learning solvers [15, 16, 17]. Learning offers a way to amortize this cost without requiring labeled optimal solutions: physics-informed networks embed governing equations directly in the training objective [18, 19], and differentiable simulation trains control policies end-to-end through the system dynamics [20]. Operator learning extends this paradigm from individual problem instances to families of them: neural operators learn mappings between functional inputs and outputs, enabling amortized inference across problem instances [21, 22, 23, 24, 25, 26]. Distribution-driven control methods have demonstrated that differentiable particle simulations can be used to train neural operators that map initial distributions to target configurations without explicitly solving MFG equations [27, 28]. Here we introduce a two-stage planner for collision-aware trajectory planning of spacecraft swarms in dense debris fields: a self-supervised, permutation-equivariant, time-conditioned neural operator, followed by a lightweight per-agent Gauss–Newton finish that closes each predicted trajectory onto exact two-body dynamics. The operator learns distribution-to-trajectory maps for spacecraft evolving under Keplerian orbital dynamics, conditions jointly on initial orbital distributions, target configurations, debris fields, and mission duration, and returns trajectories for the entire swarm through a single forward inference step. Unlike supervised operator-learning approaches that require precomputed optimal trajectories or numerical solution data, the operator is trained directly from physics-informed objectives, including terminal accuracy, fuel-effort surrogates, and closest-point-of-approach penalties, with no trajectory-optimization labels. To expose the model to informative close-approach regimes, we also generate adversarial debris by placing passive objects along nominal debris-unaware rollouts, creating training and evaluation cases in which collision avoidance requires conditioning on the debris field. This distributional viewpoint also explains why the method can generalize across swarm sizes: the operator acts on empirical samples from the underlying spacecraft, target, and debris distributions, applying the same permutation-equivariant interaction rule across arbitrary numbers of sample points. At deployment, collision-aware planning is reduced to a single forward inference pass over the N spacecraft and M debris objects, followed by a Gauss–Newton finish that is batched across the swarm; together they avoid the repeated solution of large nonlinear or mixed-integer trajectory-optimization programs as the swarm and debris field change. Our method is summarized in Fig. 1 and in Supplementary Video 1. We evaluate the method on real ephemeris data, spanning both ambient debris fields and adversarial close-approach scenarios. The planner produces low-cost, collision-aware, dynamically consistent trajectories and extrapolates stably to swarms of N=1000N=1000, providing a fast inference-time alternative to conventional optimal control in crowded orbital environments. Results We evaluate whether the learned neural operator can (i) generate low-cost orbital transfers, (i) avoid close approaches in dense debris fields, and (i) generalize beyond the training distribution in both swarm size and debris density. The operator is trained on short-duration missions of 1–12 hours, with agent counts N∈1,…,10N∈\1,…,10\ and one adversarially placed debris object per agent. We report both interpolation performance within this regime and extrapolation to swarms of up to N=1000N=1000 agents amid the full >11,000>11,000-object catalog. We evaluate the planner over a 2×22× 2 family of test conditions in which every scenario draws its N spacecraft initial states from the real Two-Line Element (TLE) catalog, so no reported result depends on synthetic starting geometry. The maneuver axis sets the retargeting magnitude: a minor maneuver perturbs each orbital element by exactly 1%1\% of its scale (random sign per element) and represents station-keeping or fine retargeting, while a major maneuver perturbs each element by exactly 10%10\% and represents rapid-response or defensive retargeting; fixing the maneuver magnitude makes fuel costs directly comparable across trials within a class. The debris axis sets the threat construction: a debris scenario surrounds the swarm with ambient catalog objects, the routine collision-avoidance regime, whereas an adversarial scenario places one worst-case object on each method’s own debris-unaware predicted path, so a planner that does not condition on the debris field is struck unless it actively deviates. These two axes span the method’s intended dual use: passive debris avoidance under realistic catalog density, and worst-case defensive pursuit-evasion against a threat that targets the nominal trajectory. Proximity is reported as a per-spacecraft rate, the percentage of planned maneuvers that pass within 100100 m of another agent or debris object, and all reported performance values are medians over 500500 Monte Carlo trials per cell. Adversarial and debris-field trajectory planning Table 1 reports terminal accuracy, fuel cost, and per-spacecraft proximity across the four scenarios and swarm sizes. We compare the operator-warm Gauss–Newton finish (GNw), which closes each operator rollout onto exact two-body dynamics, with the same finish cold-started without the operator seed (GNc). GNc reaches the same target orbit but lacks the operator’s learned collision avoidance, so the GNw–GNc gap isolates what the operator contributes. Both run batched across the swarm at every size, including N=1000N=1000, where a per-agent nonlinear program is intractable. The debris scenarios are the realistic operating regime: real spacecraft initial conditions together with the actual catalogued debris field, the setting a swarm would face in today’s crowded low Earth orbit. Here the operator-warm finish keeps close approaches rare, at or below 0.28%0.28\% of maneuvers within 100100 m, while the debris-blind GNc approaches within 100100 m on 22–3.5%3.5\% of maneuvers, an order of magnitude more, at a fuel cost within ten percent of GNw’s. This behavior holds far beyond the N≤10N≤ 10 training range, degrading gradually rather than abruptly out to N=1000N=1000 controlled spacecraft planned amid the full >11,000>11,000-object catalog, with terminal error held at 10−410^-4–10−2%10^-2\% throughout. The same learned avoidance extends to an adversarial setting. An adversarial object is one positioned directly on a spacecraft’s intended path, so that failing to deviate means near-certain collision, whether a fragment on an unusually high-conjunction trajectory or a hostile satellite maneuvering to intercept in a contested scenario. We construct this worst case by seeding one such object on each method’s own debris-unaware path. A planner blind to it is struck almost every time (GNc, 99.399.3–99.8%99.8\% of maneuvers at all sizes), whereas the operator-warm finish clears it on essentially every maneuver (at most 0.21%0.21\% at N=1000N=1000, essentially none of which involves the threat itself) at comparable fuel and accuracy. Nominal debris avoidance and defensive evasion are therefore one capability, driven by the same conditioning on the surrounding object field and exercised here against a deliberately harder threat. Table 1: Terminal accuracy, fuel cost, and collision behavior across debris and adversarial environments Error (%) Δv v (km/s) Prox100 (%) Scenario N GNw GNc GNw GNc GNw GNc Debris, minor 1 0.0007 0.0007 0.585 0.581 0 2.60 10 0.0007 0.0007 0.590 0.583 0 3.48 100 0.0007 0.0006 0.590 0.582 0.09 3.45 1000 0.0007 0.0006 0.589 0.582 0.28 3.50 Debris, major 1 0.0186 0.0140 6.004 5.592 0 2.00 10 0.0206 0.0148 6.061 5.733 0 3.42 100 0.0197 0.0148 6.084 5.826 0.05 3.35 1000 0.0197 0.0146 6.104 5.815 0.11 3.30 Adversarial, minor 1 0.0007 0.0007 0.592 0.583 0 99.8 10 0.0007 0.0007 0.589 0.582 0 99.5 100 0.0007 0.0006 0.589 0.582 0.03 99.6 1000 0.0007 0.0006 0.589 0.582 0.21 99.6 Adversarial, major 1 0.0210 0.0176 6.101 5.633 0 99.6 10 0.0213 0.0148 6.060 5.735 0 99.3 100 0.0197 0.0139 6.082 5.827 0.01 99.5 1000 0.0199 0.0151 6.090 5.873 0.07 99.5 In the single-agent setting (N=1N=1) we additionally solve the full optimal-control problem with IPOPT[29] as a reference: a nonlinear program over the controlled two-body dynamics that minimizes fuel while driving the agent to its target orbit and enforcing hard debris-avoidance constraints. It is tractable only when both the agent and debris counts are small, so we report it at N=1N=1 in the adversarial scenarios, where each agent faces a single worst-case object; against the full catalog it is intractable even at N=1N=1, as it imposes a separate collision constraint at every debris object and time step. On the real single-agent transfers it attains 0.0012%0.0012\% terminal error at Δv=0.582 v=0.582 km/s for the minor maneuver and 0.045%0.045\% at 5.525.52 km/s for the major maneuver. The operator-warm finish GNw approaches this optimum, reaching comparable terminal accuracy (Table 1) at a fuel cost within 2%2\% on the minor maneuver and within 11%11\% on the major maneuver, while remaining batched and scalable to the swarm sizes where IPOPT cannot run. IPOPT serves as a quality benchmark rather than a scalable baseline: a per-agent nonlinear program becomes prohibitively expensive as N grows, whereas the operator is evaluated once per scenario for the entire swarm. This asymmetry is reflected in Table 1: although the operator is trained only on N≤10N≤ 10, the finished trajectories retain bounded terminal errors and low proximity-violation rates at N=100N=100 and N=1000N=1000, where direct nonlinear optimization is no longer practical as a routine online planner. Figure 2 compares learned single-agent transfers with nonlinear optimal-control solutions and illustrates the effect of debris conditioning in an adversarial setting. The learned trajectories recover the optimized transfer geometry for both minor and major maneuvers. When adversarial debris is provided to the operator, the predicted trajectory changes to preserve separation above the 100100 m threshold; when the same debris information is omitted, multiple close approaches are predicted. The debris input thus directly shapes the trajectory toward avoidance. Figure 2: Learned transfer geometry and debris-conditioned avoidance. Top: Learned single-agent transfers closely match nonlinear optimal-control solutions for minor and major maneuvers. Middle: Under identical initial conditions, omitting debris information produces predicted close approaches, whereas debris conditioning keeps separation above the 100100 m safety threshold. Bottom: Log-scaled relative-trajectory spheres show the same effect geometrically: the nominal rollout intersects the debris-centered danger region, while debris conditioning reshapes the trajectory to clear it. Although the operator is trained only on N≤10N≤ 10, the finished terminal error remains tightly bounded when extrapolating to N=100N=100 and N=1000N=1000 across all four scenarios (Table 1), indicating stable degradation rather than abrupt failure far beyond the training distribution. All scenarios draw their initial conditions, and the debris scenarios their debris fields, from a publicly available Two-Line Element (TLE) catalog downloaded from Space-Track (https://w.space-track.org) on March 24, 2025, which contains over 11,00011,000 resident space objects. Each TLE is propagated to its epoch using the SGP4 orbit propagation model [30] and converted to Cartesian state vectors. The finished median terminal error stays in the 10−310^-3–10−2%10^-2\% range across swarm sizes spanning three orders of magnitude, and the N=1000N=1000 debris case plans 1,0001,000 controlled spacecraft amid the full catalog of more than 11,00011,000 objects, far outside the training distribution; collision behavior is reported in Table 1. Dynamic feasibility via a Gauss–Newton finish The operator outputs element-space trajectories that are not, by construction, hard-constrained to be the integral of a physical control sequence under exact two-body dynamics. We close each rollout with a per-agent Gauss–Newton finish: a single-shooting step over the control sequence, so every iterate is dynamically exact by construction, whose residual matches the orbit’s conserved angular-momentum and eccentricity vectors, a phase-free, non-singular terminal target, under a fuel regularizer. Warm-started from the operator (GNw) the finish inherits the operator’s collision-aware geometry; cold-started (GNc) it reaches the same target orbit without that geometry. Because the per-agent solve reduces to a small fixed-size linear system, the finish is batched across the swarm and runs at N=1000N=1000, where a per-agent nonlinear program is infeasible; its wall-clock cost is set by the trajectory length rather than by the swarm or debris count, so it stays near-constant as the swarm grows (Runtime and scaling). The full construction, including the conserved-vector residual and its batched solution, is derived in Methods. GNw drives terminal error to 10−310^-3–10−2%10^-2\% (Table 1), one to two orders below the operator’s raw element-space output and comparable to the single-agent optimum, at a fuel cost close to that single-agent optimum (Table 1), and the accuracy is bounded across three orders of magnitude in N. Collision avoidance Table 1 reports medians, in this section we report collision behavior with greater granularity for three quantities across the debris and adversarial scenarios: the operator’s raw output (ML), the operator-warm finish (GNw), and the debris-blind cold finish (GNc). The metric is the per-spacecraft rate at 100100 m. The contrast is decisive and holds at every swarm size. Against the adversarial threat the debris-blind GNc is struck on almost every maneuver (99.399.3–99.8%99.8\%), whereas the operator clears it essentially always (at most 0.21%0.21\%); under ambient catalog debris the same ordering holds at lower levels, with GNc approaching a debris object an order of magnitude more often than GNw. Throughout, GNw tracks ML to within a few tenths of a percent: the dynamics-closing finish preserves the operator’s learned avoidance rather than eroding it, and the large GNw–GNc gap is precisely the avoidance value the operator supplies through the warm start. Splitting the residual by pair type shows that avoidance of external objects is effectively complete. In the adversarial scenarios GNw’s agent–debris rate is at most 0.002%0.002\% at every swarm size (1010 events across 500,000500,000 maneuvers in the worst cell), and under catalog debris it stays at or below 0.06%0.06\% at 100100 m and 0.02%0.02\% at the tighter 5050 m radius, against 3.2%3.2\% (at 100100 m) and 12%12\% (at a conservative 500500 m screen) for GNc. What remains for ML and GNw is almost entirely agent–agent, appearing only once the swarm is dense (N≥100N≥ 100) and reflecting the swarm’s own internal density. This component is directly shaped by training: the scenario construction deliberately forces sustained agent–agent conflict (clustered starts and converging targets, Methods) under a dedicated agent–agent penalty, and the resulting operator reduces swarm-internal proximity three- to eight-fold relative to the avoidance-free GNc (0.070.07–0.22%0.22\% versus 0.510.51–0.67%0.67\% at N=1000N=1000 across the four scenarios). Further improving agent–agent deconfliction, through both training design and architecture design, is a direction for future work. (A per-trial rate, which flags an entire trial if any single agent comes close, is naturally higher and grows with swarm size, so we report the per-spacecraft rate as the true per-vehicle risk.) Figure 3 shows the distribution of each spacecraft’s closest approach for all three methods: under the adversarial scenarios GNc collapses onto the threat (∼10 10 m) while the bulk of the ML and GNw mass sits kilometers above the 100100 m threshold, and for the operator methods the only sub-threshold mass is a thin swarm-internal tail in the densest cells, consistent with the sub-0.3%0.3\% per-spacecraft rates in Table 1. Figure 3: Distribution of each spacecraft’s closest approach (500500 Monte Carlo trials per cell; one sample per spacecraft per trial). Each cell shows three violins per swarm size, for the operator output (ML, green), the operator-warm finish (GNw, red), and the cold debris-blind finish (GNc, blue); black bars mark medians and the dotted and dashed lines mark the 100100 m and 500500 m thresholds. GNc sits lowest in every cell and collapses onto the threat under the adversarial scenarios, whereas the ML and GNw mass sits well above the threshold at every swarm size; their only sub-threshold mass is a thin swarm-internal (agent–agent) tail in the densest cells, matching the per-spacecraft proximity rates in Table 1. Grid-independent proximity scoring Every separation reported in this paper is measured with a closest-point-of-approach (CPA) refinement of the sampled trajectory defined in equation (7): within each rollout interval we solve analytically for the instant of minimum separation, rather than reading off the smallest distance at the grid points. This matters because two objects usually pass closest between samples, so a grid-only scorer both overstates how far apart they get and under-counts the close approaches that fall between samples, by an amount that depends on how finely the trajectory happens to be discretized. Figure 4 makes the effect concrete. Re-scoring the densest (N=1000N=1000) cells while sweeping the grid from coarse to fine, the grid-sampled minimum overstates the true clearance by roughly 1.51.5–2×2× and drifts with resolution, so its apparent safety margin is largely an artifact of the step count. The CPA measure, by contrast, is essentially flat across the sweep, recovering the same physical closest-approach distance at every resolution. The proximity rates reported throughout are therefore both faithful, reflecting the distances objects actually reach, and grid-independent, so they cannot be made to look safer by integrating on a coarser grid. This behavior is representative: the same flatness holds across all four scenarios and every swarm size. Figure 4: CPA scoring is faithful and grid-independent. Median minimum separation for the densest debris cell (N=1000N=1000, minor maneuver) as the scoring grid is swept from coarse to fine (integration steps K; log scale). The grid-sampled minimum (orange, dashed) overstates the true clearance and drifts with resolution, only approaching the correct value as the grid is refined, whereas the within-interval CPA measure (blue, solid) used throughout is flat across the sweep, recovering the same physical closest-approach distance at every resolution. The shaded band is clearance that the grid scorer reports but that does not physically exist; the vertical line marks the grid used elsewhere in the paper. The same behavior holds across all scenarios and swarm sizes. Runtime and scaling We measured wall-clock time for both stages on the same workstation used for the Monte Carlo evaluations (a single NVIDIA RTX 2080 Ti, 11 GB). Both the operator rollout and the Gauss–Newton finish use a constant physical timestep Δt=120 t=120 s, so the number of integration steps scales with mission duration (K≈T/ΔtK≈ T/ t, from ≈31≈31 for a 1-hour transfer to ≈361≈361 for the full 12-hour horizon) rather than with swarm size. Figure 5 reports the median time per (N,M)(N,M) cell at a representative 6-hour horizon (K≈181K≈ 181). Figure 5: Runtime scaling of the two pipeline stages. Median wall-clock time (s) on a single NVIDIA RTX 2080 Ti as a function of swarm size N and debris count M. Left: neural-operator inference, which grows only mildly (0.60.6–2.92.9 s). Right: the Gauss–Newton finish, near-uniform at ≈35≈\!35 s regardless of N or M (all cells span only 3434–3636 s), since its cost is set by the per-iteration rollout length rather than the agent or debris count. Both panels use a zero-based color scale, so the finish reads as a near-flat band. Neither stage scales combinatorially. Operator inference stays below 3 seconds throughout, rising only from 0.70.7 s at N=1N=1 to 2.92.9 s at N=1000N=1000 with M=1000M=1000 debris. The Gauss–Newton finish reduces to a batched, fixed-size linear solve per agent (Methods), so its cost is set by the sequential RK4 rollouts in each iteration rather than by N or M: it runs in ≈35≈35 s and is essentially flat in swarm size (34.934.9 s at N=1N=1, 36.436.4 s at N=1000N=1000). The finish dominates, so the full pipeline replans the entire swarm in under a minute at the 6-hour horizon and roughly twice that at 12 hours; neither stage scales combinatorially with agent or debris count. The single-agent IPOPT solve, by comparison, costs ≈5.7≈5.7 s per agent and does not amortize across the swarm, serving as a quality reference at N=1N=1 rather than a scalable planner. Duration generalization The operator is trained on mission durations up to 12 hours and is not expected to extrapolate reliably to substantially longer horizons in a single rollout. The measured inference latency suggests that the operator could be embedded in a closed-loop receding-horizon controller, repeatedly querying the model with updated spacecraft and debris states while remaining within the temporal support of the training distribution. We do not evaluate this mode in the present work. All reported metrics are from single-shot rollouts on the trained horizon; long-horizon receding-horizon execution is therefore a deployment hypothesis rather than a measured result. Discussion This work demonstrated that collision-aware trajectory planning for an entire spacecraft swarm can be amortized into a single forward pass of a permutation-equivariant neural operator, trained without optimal-trajectory labels from self-supervised physics objectives and adversarial threats generated against the model’s own rollouts. Where classical planners re-solve a nonlinear program per agent and per scenario, the operator absorbs that cost at training time: trained on ten spacecraft, it transfers zero-shot to swarms two orders of magnitude larger amid the full catalogued debris field, with rare proximity violations and stable terminal accuracy, suggesting that it captures structure in the underlying dynamics rather than memorizing fixed configurations. Two design choices carry the result and are not specific to astrodynamics. The first is the division of labor between learning and numerics: the operator supplies collision-aware geometry in one batched inference, and the Gauss–Newton finish closes it onto exact dynamics, so the learned component is never asked to guarantee physics and the numerical component is never asked to see debris. The second is that the training distribution, as much as the model, determines what is learned: effective avoidance emerged from deliberately constructed scenarios, crossing-orbit threats with well-conditioned avoidance gradients and clustered, converging swarm geometries that force sustained agent–agent conflict. Several limitations remain. For very small maneuvers the finished Δv v is mildly suboptimal, because the smooth quadratic control surrogate biases the warm start away from the sharp impulse-like profiles that minimize fuel for near-zero terminal change, and accuracy degrades outside the trained 12-hour horizon, consistent with operator-learning models under distribution shift. The interaction penalties are soft and the finish carries no collision term, so the method offers no worst-case collision-avoidance guarantee. At the largest swarms the residual proximity is almost entirely agent–agent (Table 1); training on conflicting geometries under a dedicated agent–agent penalty reduces it several-fold relative to the avoidance-free baseline, and further improvement, through both training design and architecture design, is a direction for future work. Both training and evaluation assume deterministic two-body Keplerian motion, without J2J_2, atmospheric drag, solar-radiation pressure, third-body perturbations, state-estimation uncertainty or thrust execution error, effects that matter operationally because catalog uncertainty can be comparable to the 100100–500500 m thresholds considered here. Finally, extrapolation beyond the training distribution is observed empirically but not theoretically characterized. As orbits grow more congested, planning methods whose cost scales with a single batched inference rather than with the number of pairwise constraints will become a prerequisite for swarm autonomy, in routine collision avoidance and in the adversarial regime alike. The recipe demonstrated here, self-supervised physics losses, adversarial scenario generation and a certified numerical finish, is not specific to orbital mechanics, and offers a template for scalable, collision-aware multi-agent planning wherever swarms must move through contested, cluttered environments. Methods Problem setup and notation We consider a swarm of N controlled spacecraft operating in the presence of M unpowered debris objects over a fixed horizon T. Each spacecraft state is represented in Keplerian orbital elements, (t)=[aειΩων]⊤∈ℝ6, x(t)= bmatrixa& & & &ω&ν bmatrix ^6, where a is the semi-major axis, ε is the eccentricity, ι is the inclination, Ω is the right ascension of the ascending node, ω is the argument of periapsis, and ν is the true anomaly. Each debris object is represented similarly as (t)∈ℝ6 d(t) ^6. Spacecraft apply control accelerations in the Radial–Transverse–Normal (RTN) frame, (t)=[uRuSuW]⊤∈ℝ3, u(t)= bmatrixu_R&u_S&u_W bmatrix ^3, where uRu_R, uSu_S, and uWu_W denote the radial, along-track, and cross-track acceleration components, respectively. The target specification constrains only the first five orbital elements (a,ε,ι,Ω,ω)(a, , , ,ω), allowing free phasing in true anomaly ν. Orbital dynamics Spacecraft orbital element dynamics are described by the Gauss Variational Equations (GVE), a˙=2a2h(εsin(ν)uR+pruS), a= 2a^2h\! ( (ν)\,u_R+ pr\,u_S ), (1a) ε˙=1h(psin(ν)uR+((p+r)cos(ν)+rε)uS), = 1h\! (p (ν)\,u_R+ ((p+r) (ν)+r )\,u_S ), (1b) ι˙=rcos(θ)huW, = r (θ)h\,u_W, (1c) Ω˙=rsin(θ)hsin(ι)uW, = r (θ)h ( )\,u_W, (1d) ω˙=1hε(−pcos(ν)uR+(p+r)sin(ν)uS)−rsin(θ)hcot(ι)uW, ω= 1h \! (-p (ν)\,u_R+(p+r) (ν)\,u_S )- r (θ)h ( )\,u_W, (1e) ν˙=hr2+1hε(pcos(ν)uR−(p+r)sin(ν)uS), ν= hr^2+ 1h \! (p (ν)\,u_R-(p+r) (ν)\,u_S ), (1f) with p=a(1−ε2),r=p1+εcos(ν),h=μp,θ=ω+ν.p=a(1- ^2), r= p1+ (ν), h= μ\,p, θ=ω+ν. (2) For compactness, we also write the controlled dynamics as ˙=()+(), x= f( x)+ g( x) u, (3) where f captures the uncontrolled evolution and g collects the control-affine coefficients implied by equation (1). Debris are modeled as unpowered (≡ u≡ 0) and propagated under the corresponding uncontrolled dynamics ˙=() d= f( d). Distributional formulation Let P0P_0 denote the initial spacecraft distribution over ℝ6R^6, P1P_1 the target distribution over ℝ5R^5 (for the first five elements), and PdP_d the initial debris distribution over ℝ6R^6. We seek a trajectory map F(,t)F( x,t) whose pushforward defines the time-varying swarm distribution, Pt=F(⋅,t)#P0.P_t=F(·,t)_\#P_0. (4) Rather than solving a large coupled optimal control problem directly for each (P0,Pd,P1,T)(P_0,P_d,P_1,T) instance, we learn an operator that amortizes this mapping across scenarios. Neural trajectory operator We define a time-conditioned neural operator θ:(P0,Pd,P1,t,T)↦F(⋅,t)#P0,G_θ:(P_0,P_d,P_1,t,T) F(·,t)_\#P_0, (5) implemented as a permutation-equivariant transformer with multi-head attention. Our architecture builds upon the operator-learning framework introduced in Huang et al. [27], which uses sampling-invariant, permutation-equivariant attention blocks to learn solution operators for mean-field games from samples of the initial and terminal distributions. In contrast, our modified network introduces two key differences tailored to the spacecraft-swarm navigation setting: first, we include dual cross-attention streams in parallel, one between agents and the target orbits, and one between agents and the debris set, rather than a unified attention block over all inputs. This separation enables dedicated feature pathways for target alignment and debris avoidance. We maintain a shallow projection head that concatenates fused agent features, target-attention features, debris-attention features, and the original lifted queries, to output per-agent orbital-element predictions. Detail on the composition of the operator is presented in Fig. 6. Figure 6: Architecture of solution operator network, where X0X_0, X1X_1, and XdX_d are sequence of spacecraft/debris locations with arbitrary sequence length, MLP stands for Multi-layer Perceptron feed-forward network, and MHCA stands for Multi-head Cross Attention network. Terminal loss Given samples 0,ii=1N∼P0\ x_0,i\_i=1^N P_0 and associated target samples des,ii=1N∼P1\ x_des,i\_i=1^N P_1, the operator produces terminal states i(T)=F(0,i,T) x_i(T)=F( x_0,i,T). We penalize mismatch in the first five elements using (F(⋅,T)#P0;P1)=1N∑i=1N‖Π5(i(T))−des,i‖22,T\! (F(·,T)_\#P_0;P_1 )= 1N _i=1^N \| _5\! ( x_i(T) )- x_des,i \|_2^2, (6) where Π5(⋅) _5(·) projects onto (a,ε,ι,Ω,ω)(a, , , ,ω). Interaction penalties Collision avoidance is enforced using a closest-point-of-approach (CPA) penalty evaluated in Cartesian space. A state x expressed in Keplerian elements is converted to its Earth-centered inertial (ECI) position–velocity state through the standard Keplerian element-to-Cartesian map Γ , written =Γ()=(,) χ= ( x)=( p, v), and this conversion is applied to every spacecraft and debris state at every timestep. Debris are propagated as unpowered objects by advancing only the true anomaly ν under two-body Keplerian motion, keeping (a,ε,ι,Ω,ω)(a, , , ,ω) fixed. Given a spacecraft state s=(s,s) χ_s=( p_s, v_s) and a debris state d=(d,d) χ_d=( p_d, v_d) at the start of a discrete interval of duration Δt t, we compute the relative position and velocity Δ=s−d,Δ=s−d. p= p_s- p_d, v= v_s- v_d. The CPA time is tcpa(s,d)=clip(−Δ⊤Δ‖Δ‖22+ϵ, 0,Δt),t_cpa( χ_s, χ_d)=clip\! (- p v\| v\|_2^2+ε,\,0,\, t ), and the CPA separation is dcpa(s,d)=‖Δ+Δtcpa(s,d)‖2.d_cpa( χ_s, χ_d)= \| p+ v\,t_cpa( χ_s, χ_d) \|_2. (7) We penalize violations of a pair-type-dependent safety radius using a quadratic hinge, evaluated at every interval along the horizon, CPA(Ps,Pd)= _CPA(P_s,P_d)= s∼Ps,d∼Pd(κmax0,rd−dcpa(Γ(s),Γ(d)))2 _ x_s P_s,\, x_d P_d (κ \! \0,\;r_d-d_cpa\! ( ( x_s), ( x_d) ) \ )^2 +wss,s′∼Ps(κmax0,rs−dcpa(Γ(s),Γ(s′)))2, +w_s\,E_ x_s,\, x _s P_s (κ \! \0,\;r_s-d_cpa\! ( ( x_s), ( x _s) ) \ )^2, (8) where PsP_s and PdP_d are the spacecraft and debris populations, and the second expectation excludes the self-pair s′=s x _s= x_s. The spacecraft–debris safety radius rd=1r_d=1 km applies at all times; the spacecraft–spacecraft safety radius rs=100r_s=100 m applies only for t≥0.2Tt≥ 0.2\,T and its term carries the weight ws=10w_s=10; κ is a scaling constant. The spacecraft–spacecraft radius matches the 100100 m threshold used in evaluation, and the grace window over the first 20%20\% of the transfer exempts the deliberately conflicting clustered start (see Scenario Sampling), which is unavoidable by construction, while still penalizing any conflict the swarm has not resolved by mid-transfer, including the converging arrival. This yields a smooth, differentiable surrogate that concentrates gradient signal near predicted close-approach events. Adversarial debris generation. To stress-test avoidance behavior during training, every debris object is generated adversarially as a crossing-orbit threat to the model’s nominal rollout, the trajectory produced when the model is prompted with an empty debris set. Concretely, for a given (P0,P1,T)(P_0,P_1,T), we first perform a rollout with Pd=∅P_d= to obtain each agent’s nominal trajectory. For each debris object we select an agent and a random hit time thitt_hit in the latter half of the transfer, and take the agent’s nominal ECI state (,)( p, v) at thitt_hit. The debris velocity is the agent’s velocity rotated by a random angle drawn from [20∘,75∘][20 ,75 ] about a random axis perpendicular to it; the rotation preserves speed, so the debris orbit remains bound, but crosses the agent’s path with a substantial relative velocity, as in a genuine conjunction, rather than trailing it co-orbitally. The debris position is offset from the agent’s by a sub-safety-radius near-miss distance (0.350.35–0.65rd0.65\,r_d, random direction), the perturbed state is converted back to orbital elements, and, because the debris propagation model advances only ν, the initial true anomaly ν0 _0 is chosen so that the object reaches this state at thitt_hit under Keplerian motion. The near-miss offset is essential: a threat placed at exact position–velocity coincidence produces a closest-approach distance of zero, at which the CPA penalty’s gradient with respect to position vanishes identically, so the model receives a large penalty but no direction in which to evade. The crossing near-miss instead yields a well-conditioned avoidance gradient on every sample. Fuel-cost surrogate from Gauss variational dynamics To encourage fuel-efficient transfers, we penalize the magnitude of the control acceleration implied by the Gauss Variational Equations (GVE). The element-rate dynamics take the control-affine form ˙=()+() x= f( x)+ g( x)\, u, where ∈ℝ6 x ^6 is the orbital-element state, ∈ℝ3 u ^3 is the control acceleration in the rotating RTN frame, ∈ℝ6 f ^6 is the uncontrolled Keplerian rate (nonzero only in the ν component), and ()∈ℝ6×3 g( x) ^6× 3 collects the GVE control-affine coefficients. Given a predicted trajectory (t) x(t), we estimate ˙(t) x(t) by central finite differences over the discretized rollout and infer the corresponding control via the Moore–Penrose pseudoinverse † g : ^(t)=((t))†(˙(t)−((t))). u(t)= g( x(t)) ( x(t)- f( x(t)) ). We then define the fuel-cost surrogate as the time integral of squared control magnitude, ℱ=∫0T‖^(t)‖22t,F= _0^T\| u(t)\|_2^2\,dt, (9) implemented in discrete time using the rollout grid. This term does not require solving an optimal control problem during training, yet it preserves the physical meaning of minimizing RTN control effort under the element-rate dynamics. Trajectory parameterization with a biased baseline and learned residual We represent each spacecraft trajectory in orbital elements as a smooth baseline transfer augmented by a learned residual. Let τ=t/T∈[0,1]τ=t/T∈[0,1] denote normalized time, and let (t)=[a,ε,ι,Ω,ω,ν]⊤ x(t)=[a, , , ,ω,ν] . For the first five “slow” elements (a,ε,ι,Ω,ω)(a, , , ,ω), we define a deterministic baseline ¯1:5(τ) x_1:5(τ) that interpolates between the initial and target values, and we learn an additive bias (residual) Δ1:5(τ) x_1:5(τ): 1:5(τ)=¯1:5(τ)+Δ1:5(τ). x_1:5(τ)= x_1:5(τ)+ x_1:5(τ). The baseline uses linear interpolation for (a,ε,ι)(a, , ) and wrapped interpolation for angular elements (Ω,ω)( ,ω) via angdiff(α,β)=atan2(sin(α−β),cos(α−β)),angdiff(α,β)=atan2( (α-β), (α-β)), so that angle differences remain on the principal branch. The resulting slow-element sequence is projected by clamping physical bounds (e.g., a>0a>0, 0≤ε<10≤ <1, 0<ι<π0< <π) and wrapping angles to [0,2π)[0,2π). The true anomaly ν is then propagated forward in time using a physically grounded Keplerian rate with an additional learned correction. Specifically, at discrete times tkt_k with steps Δtk t_k, we update νk+1=νk+Δtk(ν˙kep(ak,εk,νk)+Δν˙k), _k+1= _k+ t_k ( ν_kep(a_k, _k, _k)+ ν_k ), where ν˙kep ν_kep is the two-body Keplerian rate and Δν˙k ν_k is the network output evaluated at the normalized time τk=tk/T∈[0,1] _k=t_k/T∈[0,1]. This parameterization biases the model toward feasible transfers while allowing the learned residual to represent nontrivial maneuver geometry and phasing behavior. Training objective For each sampled scenario (P0,Pd,P1,T)(P_0,P_d,P_1,T), the underlying planning task can be viewed as a finite-horizon collision-aware optimal control problem. Given spacecraft samples 0,ii=1N∼P0\ x_0,i\_i=1^N P_0, target samples des,ii=1N∼P1\ x_des,i\_i=1^N P_1, and debris samples 0,jj=1M∼Pd\ d_0,j\_j=1^M P_d, the ideal per-instance problem is to find trajectories and controls that minimize fuel expenditure while reaching the target distribution and avoiding close approaches: mini,ii=1N _\ x_i, u_i\_i=1^N ∑i=1N∫0T‖i(t)‖22t+λT∑i=1N‖Π5(i(T))−Π5(des,i)‖22+λICPA _i=1^N _0^T\| u_i(t)\|_2^2\,dt+ _T _i=1^N \| _5\! ( x_i(T) )- _5\! ( x_des,i ) \|_2^2+ _I\,J_CPA (10) s.t. .t. ˙i(t)=(i(t))+(i(t))i(t),i(0)=0,i, x_i(t)= f( x_i(t))+ g( x_i(t)) u_i(t), x_i(0)= x_0,i, ˙j(t)=(j(t)),j(0)=0,j, d_j(t)= f( d_j(t)), d_j(0)= d_0,j, where Π5 _5 projects onto the first five orbital elements (a,ε,ι,Ω,ω)(a, , , ,ω) and CPAJ_CPA denotes the closest-point-of-approach penalty, evaluated separately over spacecraft–debris and spacecraft–spacecraft pairs with their respective safety radii (see Interaction penalties). This formulation expresses the desired collision-aware planning problem for a single scenario, but solving equation (10) repeatedly would require expensive numerical optimization and would make supervised training dependent on a large library of precomputed optimal trajectories. Instead, we train θG_θ as an amortized solution operator over a distribution of scenarios. At each training instance, we sample (P0,Pd,P1,T)(P_0,P_d,P_1,T), construct a debris set with an adversarial fraction, and generate trajectories using the biased baseline plus learned residual parameterization, namely θ(t;0)=¯(t)+Δθ(t;0), x_θ(t; x_0)= x(t)+ x_θ(t; x_0), (11) where ¯ x is the deterministic baseline of Trajectory parameterization with a biased baseline and learned residual and Δθ x_θ is the network residual, conditioned on the full scenario (P0,Pd,P1,T)(P_0,P_d,P_1,T). The rollout realizes the trajectory map of equation (5) sample-wise: θ(P0,Pd,P1,t,T)=θ(t;⋅)#P0G_θ(P_0,P_d,P_1,t,T)= x_θ(t;·)_\#P_0. The network parameters are optimized by minimizing the expected self-supervised loss minθ(P0,Pd,P1,T)∼[λfℱθ+λIℐθ+λTθ+λSθ], _θ\;E_(P_0,P_d,P_1,T) [ _fF_θ+ _II_θ+ _TT_θ+ _SS_θ ], (12) where D is the training distribution over planning scenarios. The four terms are defined as follows: ℱθ=0∼P0∫0T‖(θ(t;0))†(˙θ(t;0)−(θ(t;0)))‖22tF_θ=E_ x_0 P_0 _0^T \| g( x_θ(t; x_0)) ( x_θ(t; x_0)- f( x_θ(t; x_0)) ) \|_2^2\,dt (13) is the fuel cost, the distributional form of equation (9); and ℐθ=CPA(θ(t;⋅)#P0,Pd)I_θ=J_CPA ( x_θ(t;·)_\#P_0,\;P_d ) (14) is the closest-point-of-approach interaction penalty of equation (8), evaluated on the pushforward of P0P_0 under the rollout along the horizon; and θ=(0,des)∼(P0,P1)‖Π5(θ(T;0))−Π5(des)‖22T_θ=E_( x_0, x_des) (P_0,P_1) \| _5\! ( x_θ(T; x_0) )- _5\! ( x_des ) \|_2^2 (15) is the terminal loss of equation (6), where the expectation is over paired samples, each spacecraft with its own assigned target; and θ=0∼P0‖Π5(θ(0;0))−Π5(0)‖22S_θ=E_ x_0 P_0 \| _5\! ( x_θ(0; x_0) )- _5\! ( x_0 ) \|_2^2 (16) is the initial loss, which anchors the rollout at τ=0τ=0 to the sampled initial state; in both anchoring terms the projection Π5 _5 excludes the true anomaly. This objective trains the operator without ground-truth optimal trajectories or numerically generated trajectory labels, while preserving the structure of the per-instance optimal control problem in equation (10). The specific choice of the weights (λf,λI,λT,λS)( _f, _I, _T, _S) is discussed in Table 5. Dynamic-feasibility finish via Gauss–Newton terminal targeting The operator’s element-space rollout is not, by construction, the integral of a physical control sequence under exact two-body dynamics. We finish each rollout, per agent, with a single-shooting Gauss–Newton (GN) step with Levenberg–Marquardt damping. The control sequence 0:T−1 u_0:T-1 is the only decision variable; the state is the exact RK4 rollout k+1=ΦΔtRK4(k,k) x_k+1= ^RK4_ t( x_k, u_k) from the fixed initial state, so every iterate is dynamically feasible and there are no dynamics constraints. We target the orbit through its conserved vectors rather than its element angles: the specific angular momentum =× h= r× v and the eccentricity vector =(×)/μ−/‖ e=( v× h)/μ- r/\| r\|. Both are invariant along an orbit, so the residual is phase-free, and smooth in (,)( r, v), so the Jacobian stays well-defined for the near-circular and near-equatorial orbits where element-angle residuals are singular. With targets (⋆,⋆)( h , e ) from the goal orbit P1P_1, the residual and objective are ()=[((T)−⋆)/sh;((T)−⋆)/se]∈ℝ6,min‖()‖22+λf‖22. ρ( u)= [( h(T)- h )/s_h;\;( e(T)- e )/s_e ] ^6, _ u~\| ρ( u)\|_2^2+ _f\| u\|_2^2. (17) Each iteration forms the Jacobian =∂/∂∈ℝ6×3T J=∂ ρ/∂ u ^6× 3T by automatic differentiation through the RK4 rollout, six vector–Jacobian products, one per residual component, sharing a single retained backward graph, and takes the damped Gauss–Newton step Δ=−(⊤+(λf+μ))−1(⊤+λf), u=- ( J J+( _f+μ) I )^-1 ( J ρ+ _f u ), (18) where μ≥0μ≥ 0 is the Levenberg–Marquardt damping. The decision vector u has dimension 3T3T (hundreds to thousands), but the residual has only six components, so we never form the 3T×3T3T× 3T system. Writing α=λf+μα= _f+μ and =−(⊤+λf) b=-( J ρ+ _f u), the Woodbury identity collapses equation (18) to a single 6×66× 6 solve, Δ=1α(−⊤(α6+⊤)−1), u= 1α ( b- J (α I_6+ J J )^-1 J\, b ), (19) whose cost is independent of the horizon length T and identical for every agent. The swarm is therefore solved as one batched stack of 6×66× 6 systems, the property that lets the finish run at N=1000N=1000, where a per-agent nonlinear program is intractable. The 6×66× 6 inner solve is carried out in double precision to absorb the 1/α1/α cancellation when the fuel term is active (λf>0 _f>0); for the pure terminal target (λf=0 _f=0) the step reduces to the numerically benign minimum-norm form Δ=−⊤(⊤+μ6)−1 u=- J ( J J +μ I_6)^-1 ρ. The damping is adapted by a per-agent trust region. At each iteration the trial step is accepted for every agent whose combined objective ‖22+λf‖22\| ρ\|_2^2+ _f\| u\|_2^2 decreases and rejected otherwise; μ is scaled down by 0.50.5 on acceptance and up by 44 on rejection, clamped to a fixed range. Because acceptance is decided independently per agent inside the batched iteration, a stiff agent that needs heavy damping does not stall the rest of the swarm. We run at most 2525 iterations and stop early once the terminal residual falls below 10−710^-7. The same conserved-vector targeting defines the single-agent IPOPT reference. We report two variants. GNw warm-starts u from the operator rollout via an analytical RK4 inverse, per interval, inverting the constant input-to-state map B=[(Δt2/2)I3;ΔtI3]B=[( t^2/2)I_3;\, t\,I_3] on the zero-control residual, so the finish begins from the operator’s collision-aware geometry; GNc cold-starts from =0 u=0. The fuel regularizer λf _f (a small default; λf=0 _f=0 recovers the pure minimum-norm step) trades a slight amount of terminal accuracy for lower control energy, and is added as the residual rows in equation (18) so the solve stays 6×66× 6. The finish carries no collision term: avoidance is supplied entirely by the collision-aware operator seed, and the minimum-norm character of the step perturbs that trajectory only as much as exact dynamics require. All reported Δv v and terminal-error values are computed on the finished, dynamically-exact trajectory. Architecture and training details This subsection gives the exact layer structure of the operator θG_θ (equation (5)), its initialization, and the training hyperparameters. Input feature lifting Each spacecraft or debris state =[a,ε,ι,Ω,ω,ν]⊤∈ℝ6 x=[a,\, ,\, ,\, ,\,ω,\,ν] ^6 is expressed in non-dimensional units using the scale factors Lc=6.37×106m,Vc=μE/Lc,Tc=Lc/Vc≈803s.L_c=6.37× 10^6\ m, V_c= _E/L_c, T_c=L_c/V_c≈ 803\ s. (20) Before encoding, every point in P0P_0, P1P_1, and PdP_d is lifted from ℝ6R^6 to ℝ8R^8 by appending the normalized time τ∈[0,1]τ∈[0,1] and the non-dimensional mission duration Tnd=T/TcT_nd=T/T_c: ~(j)=[(j),τ,Tnd]∈ℝ8. x^(j)= [ x^(j),\;τ,\;T_nd ] ^8. (21) This conditioning is applied identically across all three input sets, yielding batched tensors of shape (B,N,8)(B,N,8), (B,N,8)(B,N,8), and (B,M,8)(B,M,8) for spacecraft, target, and debris, respectively. Network architecture The operator is structured in three sequential stages. Per-Point Encoders Three independent two-layer pointwise MLPs, implemented as 1×11×1 convolutions over the particle dimension to preserve permutation equivariance,project each input stream into a shared latent space of width h: ϕ:ℝ8→ℝh,ϕ(⋅)=σ(W2σ(W1⋅+b1)+b2),φ:\;R^8 ^h, φ(·)=σ (W_2\,σ(W_1\,·+b_1)+b_2 ), (22) where σ denotes GELU and dropout (rate p) is applied after the first linear map. Separate parameter sets ϕ0 _0, ϕ1 _1, ϕd _d are maintained for the spacecraft, target, and debris streams, producing encoded representations 0∈ℝB×h×N,1∈ℝB×h×N,d∈ℝB×h×M.H_0 ^B× h× N, _1 ^B× h× N, _d ^B× h× M. (23) Dual Cross-Attention Stack The operator employs two parallel cross-attention streams. This appendix gives the precise per-layer update. At each layer ℓ=1,…,L =1,…,L, the spacecraft encodings q(ℓ−1)∈ℝN×B×hX_q^( -1) ^N× B× h serve as queries into both streams: 1(ℓ) _1^( ) =MHA1(q(ℓ−1),1,1), =MHA_1\! (X_q^( -1),\,H_1,\,H_1 ), (24) d(ℓ) _d^( ) =MHAd(q(ℓ−1),d,d), =MHA_d\! (X_q^( -1),\,H_d,\,H_d ), (25) q(ℓ) _q^( ) =q(ℓ−1)+1(ℓ)+d(ℓ), =X_q^( -1)+Y_1^( )+Y_d^( ), (26) where each MHAMHA block is a standard scaled dot-product multi-head attention with hheadsh_heads heads, head dimension h/hheadsh/h_heads, and no projection biases. The two streams share no parameters. The residual accumulation in the last line propagates both goal-directed and avoidance signals simultaneously through the depth of the stack. Final Decoder After L attention layers, features from the accumulated query q(L)X_q^(L), the final target-attention output 1(L)Y_1^(L), the final debris-attention output d(L)Y_d^(L), and the raw lifted inputs are concatenated per point: =[q(L),1(L),d(L),P~0,P~1]∈ℝB×(3h+16)×N,Z= [X_q^(L),\;Y_1^(L),\;Y_d^(L),\; P_0,\; P_1 ] ^B×(3h+16)× N, (27) where P~0 P_0 and P~1 P_1 are the lifted (dimension-8) inputs, and 3h+2×8=3h+163h+2×8=3h+16. A final two-layer pointwise MLP projects this to the output: Conv1d(3h+16→h)→Dropout(p)→GELU→Conv1d(h→7).Conv1d(3h+16→ h) (p) 1d(h→ 7). (28) The seven output channels per agent correspond to residuals on the five slow elements (a,ε,ι,Ω,ω)(a, , , ,ω) and one residual correction to the true-anomaly rate Δν˙ ν, consistent with the trajectory parameterization. Hyperparameters and Parameter Count Table 2 lists the architectural hyperparameters. With these settings, the approximate parameter counts by component are given in Table 3. Table 2: Architectural hyperparameters of the neural operator. Hyperparameter Value Hidden dimension h 1024 Cross-attention layers L 5 Attention heads hheadsh_heads 8 Head dimension h/hheadsh/h_heads 128 Encoder MLP depth 2 layers (per stream) Decoder MLP depth 2 layers Dropout rate p 0.05 Activation σ GELU Attention projection bias None Table 3: Approximate parameter counts by component (h=1024h=1024, L=5L=5). Component Parameters Three encoder MLPs (ϕ0,ϕ1,ϕd _0, _1, _d) ≈ 3.18×106≈\,3.18× 10^6 Ten cross-attention modules (2 streams × L layers) ≈ 41.9×106≈\,41.9× 10^6 Final decoder MLP ≈ 3.17×106≈\,3.17× 10^6 Total ≈ 48.3×106≈\,48.3× 10^6 Weight Initialization All 1×11×1 convolutional layers are initialized with Kaiming uniform initialization (a=5a= 5, fan-in mode); any standalone linear layers use Xavier uniform. All biases are initialized to zero. The final output convolution is not rescaled, so initial outputs are of unit magnitude. Inference rollout At inference, the operator is queried on a uniform grid of K+1K+1 time points, τk=k/K _k=k/K, with K determined by the physical time step Δtphys=120s t_phys=120\ s and the mission duration T. The slow-element sequence is produced in a single batched forward pass by stacking all K+1K+1 time queries, and the true anomaly ν is then integrated sequentially using the recurrence in equation (1) augmented by the network’s residual correction Δν˙k ν_k. Training hyperparameters Scenario Sampling At each training step a single LEO reconfiguration scenario is drawn. The number of controlled spacecraft N∈1,…,10N∈\1,…,10\ is sampled uniformly at random. Initial states are generated as a dense cluster: a cluster-center orbit is drawn uniformly within the bounds given in Table 4 (subject to the feasibility constraint that the periapsis exceeds Earth ++ 100 km and the apoapsis lies below Earth ++ 2000 km), converted to a Cartesian state, and each agent is placed uniformly within a 5050 m ball of the center in both position and a matched velocity spread, so pairwise separations span roughly 0–100100 m. This manufactures genuine, unavoidable agent–agent conflict at every N≥2N≥ 2; LEO coverage comes from randomizing the cluster center across scenarios rather than from scattering agents, under which agent–agent conflict essentially never arises at N≤10N≤ 10. Table 4: Orbital element sampling bounds used during training. Element Lower bound Upper bound a (semi-major axis) 6848 km 8748 km ε (eccentricity) 0.001 0.1 ι (inclination) 0 π Ω (RAAN) 0 2π2π ω (argument of periapsis) 0 2π2π ν (true anomaly) 0 2π2π Targets are per-agent but converging: one maneuver class is drawn per scenario (minor or major, equal probability), one deviation vector is sampled at that class’s magnitude, each of the five slow elements perturbed by exactly its class scale with an independent random sign and the terminal true anomaly left free, and the same deviation is applied to every agent’s own P0P_0. Each agent therefore flies an exactly-1%1\% or exactly-10%10\% maneuver of its own orbit, while the targets form a cluster mirroring the start cluster: an unmodified (non-deconflicting) rollout keeps the phase-locked formation in conflict all the way to a co-located arrival, and the escape, a phase-staggered arrival, is free because the target constrains only the five slow elements. The mission duration T is drawn uniformly from [1hr,12hr][1\ hr,12\ hr]. Debris consist of one adversarially generated crossing-orbit object per agent (M=NM=N), using the construction described in Interaction penalties. Loss Weights The complete weighted objective is given by equation (12), with coefficients given in Table 5. LEO feasibility is maintained by construction rather than by penalty: scenario sampling reflects any target perturbation that would leave the feasible set back toward its interior, and numerical guards (see Numerical Stability) protect against ill-conditioned element configurations. The terminal and initial losses T and S are normalized element-wise by the scale vector [1.37, 0.1,π, 2π, 2π][1.37,\;0.1,\;π,\;2π,\;2π] for (a,ε,ι,Ω,ω)(a, , , ,ω) to account for their differing numerical ranges. Table 5: Loss term weights used during training. Term Weight Role ℱF (fuel) λf=10−2 _f=10^-2 GVE-implied control energy. ℐI (collision) λI=103 _I=10^3 Split CPA quadratic hinge penalty, κ=104κ=10^4: spacecraft–debris with rd=1kmr_d=1\ km at all times; spacecraft–spacecraft with rs=100mr_s=100\ m, weight ws=10w_s=10, active for t≥0.2Tt≥ 0.2\,T. T (terminal) λT=102 _T=10^2 Scaled MSE at τ=1τ=1 vs. P1P_1. S (initial) λS=102 _S=10^2 Scaled MSE at τ=0τ=0 vs. P0P_0 (excluding ν). Optimizer The network is trained with Adam using the hyperparameters in Table 6. Gradients are clipped by global ℓ2 _2 norm before each parameter update. Table 6: Optimizer and scheduler hyperparameters. Hyperparameter Value Optimizer Adam Learning rate 2×10−52× 10^-5 LR schedule Reduce-on-plateau (patience 1000) Gradient clip (ℓ2 _2 norm) 100 Training iterations 14,000 Random seed 10 Floating-point precision float32 Numerical Stability Three mechanisms guard against training instabilities that arise from stiff orbital dynamics and adversarial debris generation. 1. EMA spike guard. An exponential moving average of the loss (α=0.99α=0.99) is maintained throughout training. Any step whose loss exceeds max(103ℒ¯EMA, 102) (10^3\, L_EMA,\;10^2) is skipped, and the learning rate is reduced by a factor of 0.90.9 before the next step. 2. Non-finite clamping. All network outputs and intermediate Cartesian coordinate conversions (orbital elements to inertial position–velocity for the CPA computation) are passed through nan_to_num with zero replacement to prevent gradient corruption from ill-conditioned element configurations. 3. Per-term gradient validation. Each loss term is back-propagated individually with graph retention before the full backward pass, and parameter gradients are checked for non-finite values. An infinite or NaN gradient raises a runtime exception and halts training, enabling targeted diagnosis rather than silent divergence. Declarations Funding This research received no external funding. Conflict of interest / Competing interests The authors declare no competing interests. Ethics approval and consent to participate Not applicable; this computational study involves no human participants, animals, or biological materials. Consent for publication Not applicable. Data availability The Two-Line Element catalog used for initial conditions and debris fields is publicly available from Space-Track (https://w.space-track.org; snapshot of 24 March 2025). The processed ephemeris file, trained model weights, and the Monte Carlo evaluation outputs supporting all tables and figures will be deposited in a public repository with a DOI upon publication, and are available to editors and reviewers on request during assessment. Materials availability Not applicable. Code availability The training and evaluation code used to generate all results in this study will be released in a DOI-minting public repository upon publication, and is likewise available to editors and reviewers on request. Use of artificial intelligence tools Large language model tools (Anthropic Claude) assisted with code development and manuscript editing; all methods, results, analyses, and conclusions were developed and verified by the authors. No generative artificial intelligence was used to create any figure or image in this manuscript. Author contributions S.D.S. conceived the study, developed the methodology, implemented the experiments, performed the analysis, and wrote the manuscript. R.L. and S.M. supervised the research, contributed key ideas, and revised the manuscript. Z.L. and S.G. advised on technical aspects of implementation and contributed to manuscript drafting and revision. All authors reviewed and approved the final manuscript. References European Space Agency Space Debris Office [2026] European Space Agency Space Debris Office: ESA’s Annual Space Environment Report. Technical report, ESA/ESOC, Darmstadt, Germany (2026) Space.com [2026] Space.com: Every SpaceX Starlink satellite has to dodge a collision almost weekly, and experts fear the worst. Space.com. Reporting SpaceX’s semi-annual orbital-safety filing to the FCC, covering December 2025–May 2026 (2026) Lee and Ho [2023] Lee, H.W., Ho, K.: Regional constellation reconfiguration problem: Integer linear programming formulation and Lagrangian heuristic method. Journal of Spacecraft and Rockets 60(6), 1828–1845 (2023) Basu et al. [2023] Basu, H., Pedari, Y., Almassalkhi, M., Ossareh, H.R.: Computationally efficient collision-free trajectory planning of satellite swarms under unmodeled orbital perturbations. Journal of Guidance, Control, and Dynamics 46(8), 1548–1563 (2023) Eren et al. [2017] Eren, U., Prach, A., Koçer, B.B., Raković, S.V., Kayacan, E., Açıkmeşe, B.: Model predictive control in aerospace systems: Current state and opportunities. Journal of Guidance, Control, and Dynamics 40(7), 1541–1566 (2017) Chen et al. [2024] Chen, R., Dong, M., Bai, Y., Zhao, Y., Chen, X.: Trajectory planning and control of spacecraft avoiding dynamic debris swarm. Aerospace Science and Technology 151, 109273 (2024) van den Berg et al. [2011] van den Berg, J., Guy, S.J., Lin, M., Manocha, D.: Reciprocal n-body collision avoidance. In: Robotics Research. Springer Tracts in Advanced Robotics, vol. 70, p. 3–19. Springer, Berlin, Heidelberg (2011). Pedari et al. [2023] Pedari, Y., Basu, H., Ossareh, H.R.: A novel framework for trajectory planning and safe navigation of satellite swarms. IFAC-PapersOnLine 56(2), 547–552 (2023) Jung and Chung [2025] Jung, I., Chung, D.: Genetic algorithm-based approach for improving temporal resolution in constellation operation of national satellites. International Journal of Aeronautical and Space Sciences 26(1), 314–326 (2025) Xu et al. [2024] Xu, L., Zhang, G., Qiu, S., Cao, X.: Reinforcement learning-based multi-impulse rendezvous approach for satellite constellation reconfiguration. Acta Astronautica 224, 325–337 (2024) Kuhl et al. [2025] Kuhl, W., Wang, J., Eddy, D., Kochenderfer, M.J.: Markov decision processes for satellite maneuver planning and collision avoidance. In: 2025 IEEE Aerospace Conference, p. 1–9. IEEE, Big Sky, MT (2025). An et al. [2026] An, X., Luo, S., Zhang, H., Yang, Q., Ma, Y., Wang, B., Du, J., Wang, Q.: Autonomous navigation of intelligent microrobotic swarms in unknown environments. Nature Machine Intelligence 8(6), 955–968 (2026) Bensoussan et al. [2013] Bensoussan, A., Frehse, J., Yam, P.: Mean Field Games and Mean Field Type Control Theory. SpringerBriefs in Mathematics. Springer, New York (2013). Wang et al. [2022] Wang, G., Yao, W., Zhang, X., Li, Z.: A mean-field game control for large-scale swarm formation flight in dense environments. Sensors 22(14), 5437 (2022) Guo et al. [2019] Guo, X., Hu, A., Xu, R., Zhang, J.: Learning mean-field games. In: Advances in Neural Information Processing Systems, vol. 32. Curran Associates, Inc., Red Hook, NY (2019). Ruthotto et al. [2020] Ruthotto, L., Osher, S.J., Li, W., Nurbekyan, L., Fung, S.W.: A machine learning framework for solving high-dimensional mean field game and mean field control problems. Proceedings of the National Academy of Sciences 117(17), 9183–9193 (2020) Laurière et al. [2022] Laurière, M., Perrin, S., Pérolat, J., Girgin, S., Muller, P., Élie, R., Geist, M., Pietquin, O.: Learning in Mean Field Games: A Survey. Preprint at https://arxiv.org/abs/2205.12944 (2022) Raissi et al. [2019] Raissi, M., Perdikaris, P., Karniadakis, G.E.: Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, 686–707 (2019) Karniadakis et al. [2021] Karniadakis, G.E., Kevrekidis, I.G., Lu, L., Perdikaris, P., Wang, S., Yang, L.: Physics-informed machine learning. Nature Reviews Physics 3(6), 422–440 (2021) Zhang et al. [2025] Zhang, Y., Hu, Y., Song, Y., Zou, D., Lin, W.: Learning vision-based agile flight via differentiable physics. Nature Machine Intelligence 7(6), 954–966 (2025) Lu et al. [2021] Lu, L., Jin, P., Pang, G., Zhang, Z., Karniadakis, G.E.: Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence 3(3), 218–229 (2021) Berner et al. [2026] Berner, J., Liu-Schiaffini, M., Kossaifi, J., Duruisseaux, V., Bonev, B., Azizzadenesheli, K., Anandkumar, A.: Principled approaches for extending neural architectures to function spaces for operator learning. Nature Machine Intelligence 8, 1173–1181 (2026) Wang et al. [2021] Wang, S., Wang, H., Perdikaris, P.: Learning the solution operator of parametric partial differential equations with physics-informed deeponets. Science Advances 7(40), 8605 (2021) Yang et al. [2023] Yang, L., Liu, S., Meng, T., Osher, S.J.: In-context operator learning with data prompts for differential equation problems. Proceedings of the National Academy of Sciences 120(39), 2310142120 (2023) Xiao et al. [2025] Xiao, P., Zheng, M., Jiao, A., Yang, X., Lu, L.: Quantum DeepONet: Neural operators accelerated by quantum computing. Quantum 9, 1761 (2025) Xu et al. [2025] Xu, W., Han, J., Lai, R.: Self-Supervised Amortized Neural Operators for Optimal Control: Scaling Laws and Applications. Preprint at https://arxiv.org/abs/2512.24897 (2025) Huang and Lai [2025] Huang, H., Lai, R.: Unsupervised solution operator learning for mean-field games. Journal of Computational Physics 537, 114057 (2025) Cole et al. [2026] Cole, F., Wang, D., Chen, Y., Lu, Y., Lai, R.: In-Context Operator Learning on the Space of Probability Measures. Preprint at https://arxiv.org/abs/2601.09979 (2026) Wächter and Biegler [2006] Wächter, A., Biegler, L.T.: On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming 106(1), 25–57 (2006) Vallado et al. [2006] Vallado, D.A., Crawford, P., Hujsak, R., Kelso, T.S.: Revisiting spacetrack report #3. In: AIAA/AAS Astrodynamics Specialist Conference and Exhibit. American Institute of Aeronautics and Astronautics, Keystone, CO (2006).