Multi-Timescale Context-Aware Radio Resource Optimization for Harmonious NTN-TN Coexistence: A Digital Twin and Open RAN Enabled Framework

1

Multi-Timescale Context-Aware Radio Resource Optimization for Harmonious NTN-TN Coexistence: A Digital Twin and Open RAN Enabled Framework

arXiv:2610.32827v1 [cs.NI] 23 Oct 2026

Ecem Izde Kaya, Md. Zoheb Hassan, Member, IEEE, and Georges Kaddoum, Senior Member, IEEE

Abstract—Coexistence between non-terrestrial networks (NTNs) and terrestrial networks (TNs) will be essential for future wireless systems. While spectrum sharing with TNs will enable NTNs to more efficiently expand coverage and use resources, a key challenge that remains is interference. In particular, terrestrial downlink transmissions can generate harmful aggregate interference towards low-Earth orbit (LEO) satellite receivers, making interference-aware resource management essential for coexistence. To address this need, in this article, we present DT-STABLE, an Open Radio Access Network (O-RAN)-enabled, digital-twin (DT)-assisted framework for context-aware NTN–TN coexistence. By using site-specific contextual information, the DT estimates interference perceived at the LEO satellites, enabling proactive interference management without extensive real-time measurement overhead. We formulate a joint beam activation and user association problem accounting for RF-chain limits, load balancing, and aggregate satellite interference constraints, prove its NP-hardness, and develop an efficient DT-guided algorithm for practical scheduling. DT-STABLE executes a multitimescale optimization by coupling short-timescale radio resource adaptation with periodic mid- and long-timescale updates of the DT’s predictive state. This synchronization is maintained through a continuous two-way exchange of context, measurements, and control information over standardized O-RAN interfaces, ensuring that the twin remains accurate and aligned with the physical network. The results of site-specific simulations under worst-case satellite geometry at 28 GHz band, with LEO satellites at low elevation angles, demonstrate that DT-STABLE outperforms existing approaches in balancing rate and interference utility, achieving 89.63% interference-to-noise ratio (INR) satisfaction and a mean sum rate of 210.96 bps/Hz, while attaining nearly 100% INR satisfaction when satellite protection is prioritized. The proposed framework also reduces beam measurement overhead by 99.71% as compared with exhaustive sweeping and achieves up to 33.1× lower runtime than the Gurobi-based exact optimization benchmark.

I. I NTRODUCTION The integration of non-terrestrial networks (NTNs) with terrestrial networks (TNs) has been recognized by the International Telecommunication Union (ITU) and the Third Generation Partnership Project (3GPP) as a fundamental component of future 5G and 6G ecosystems [1]. However, while this integration has a significant potential in terms of enhancing global connectivity, spectrum scarcity poses a critical challenge. E. I. Kaya and G. Kaddoum are with the École de technologie supérieure (ETS), Université du Québec, Montréal, QC, Canada. M. Z. Hassan is with the Université Laval, Québec City, QC, Canada. (Corresponding author’s e-mail: md-zoheb.hassan@gel.ulaval.ca). This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible.

As conventional sub-6 GHz bands become increasingly congested [2], cellular systems have extended into frequency range 2 (FR2), spanning 24.25–52.6 GHz [3]. At the same time, portions of this spectrum overlap with Ka-band frequencies used by satellite communication systems. In particular, the Federal Communications Commission (FCC) has authorized SpaceX’s Gen2 Starlink non-geostationary orbit (NGSO) constellation for Earth-to-space operations in the 27.5–29.1 GHz and 29.5–30.0 GHz bands [4]. This spectral overlap creates a coexistence challenge, as terrestrial downlink transmissions may cause harmful interference to satellite uplinks and degrade satellite performance [5]. While 5G base stations (BSs) are primarily designed to serve terrestrial users, practical antenna patterns can radiate part of the transmitted power towards the sky via side lobes, highlighting the need for coordinated spectrum sharing solutions. Narrow and highly directional beams at millimetre-wave (mmWave) frequencies can improve spatial isolation by steering energy toward the intended receiver and limiting emissions in other directions [6]. However, interference remains a critical challenge for the reasons outlined below. • Even with narrow beams, there are sidelobes, and sidelobe suppression broadens the main lobe [7], [8]. • Satellite interference can be significantly affected not only by the direct emissions of terrestrial transmitters, but also by the upscattering of terrestrial signals by surrounding objects such as buildings, trees, and vehicles, which further increases the aggregate interference observed at the satellite receiver [9]. • Aggregate interference remains a significant concern, as true coexistence is determined by the overall antenna pattern and network deployment strategy, rather than by the interference characteristics of a single beam and single BS. In recent years, the integration of open radio access network (O-RAN) architectures with advanced resource-management methods has gained increasing attention due to its ability to support flexible, scalable, and cost-efficient network operation [10]. O-RAN naturally supports these functions through programmable control loops, standardized interfaces, and AI/MLoriented design, all of which enable virtual network intelligence to interact with physical RAN operation [11]. Digital twin networks (DTNs) are virtual real-time replicas of physical networks supporting network design, diagnosis, simulation, what-if analysis, and AI/ML-driven optimization and control [12].The integration of digital twins (DTs) with O-RAN can improve spectrum management, enhance resource allocation,

2

and support control-policy validation in a virtual environment before physical-network deployment [13]. In this paper, we propose DT-STABLE, an O-RANenabled DT-assisted multi-time-scale radio resource management (RRM) framework for dynamic TN–NTN spectrum coexistence. Unlike conventional beam management strategies that rely on frequent real-time measurements, DT-STABLE uses a site-specific channel knowledge map (CKM), which provides location-specific channel knowledge, reflecting the local wireless propagation environment [26]. The CKM is generated offline from 3D environmental geometry, material information, and RF propagation modeling of a highfidelity DT to predict terrestrial link quality and satellite interference exposure. This enables near real-time (near-RT), proactive beam activation and user scheduling under satellite visibility, aggregate interference, RF-chain, and beam-load constraints. DT-STABLE employs O-RAN closed-loop control to detect CKM prediction mismatches in dynamic wireless environments and to selectively trigger low-overhead beam sweeping for DT refinement over multi-time scale horizon. Major contributions can be summarized as follows: • We formulate the NTN–TN spectrum coexistence problem as a constrained joint beam-activation and user-association problem that balances terrestrial service utility and satellite protection through a soft aggregate-interference penalty. We further show that the resulting problem is NP-hard. • We develop DT-STABLE, an O-RAN-enabled Digital Twin Satellite-Threshold-Aware Beam–Load Engine, where DTin-the-loop optimization is built on top of the O-RAN architecture to support multi-timescale interference-aware terrestrial network control. The proposed framework jointly integrates long-term site-specific CKM construction, midterm CKM refinement, and short-term beam–user scheduling across multiple control time scales. • By using a multi-domain context-aware DT (CDT), we develop a 3D geo-environment and material-aware CKM that captures both terrestrial downlink channels and satellite interference, and facilitates proactive beam selection without exhaustive measurements. The terrestrial CKM is constructed using high-fidelity site-specific ray tracing (RT) with NVIDIA Sionna RT [27], incorporating real-world 3D building models, actual base station (BS) locations, and beamforming gains. The non-terrestrial CKM uses the OpenNTN channel model [28], incorporating two-line element (TLE) data-based satellite positions, satellite and BS antenna gains, stochastic channel variations, and propagation losses due to shadowing, clutter, weather, and atmospheric conditions. • We introduce a low-overhead, periodic, and adaptive refinement mechanism that uses limited beam sweeping to update outdated or mismatched terrestrial CKM predictions, thereby enhancing robustness against DT modeling uncertainty. • Using comprehensive site-specific simulations over the 28 GHz band, we showcase DT-STABLE’s efficacy for practical NTN-TN (uplink-downlink) coexistence scenarios while covering site-specific propagation, satellite visibility, and benchmark comparisons; we also evaluate its robustness under imperfect DT knowledge.

The remainder of this article is organized as follows. A summary of related works is provided in Section II. Section III presents an overview of the system considered. In Section IV, we present the DT-STABLE framework. The formulation of the optimization problem, the algorithm, and the multitime scale resource optimization framework are described in Section V. Finally, the performance evaluation and conclusion are presented in Sections VI and VII, respectively. II. R ELATED W ORKS Resource Management for NTN-TN: Unlike cellular networks, NTN systems face dynamic channel and topology conditions from satellite mobility, long propagation distances, and time-varying visibility, which motivates dedicated resourcemanagement strategies. In low-Earth orbit (LEO) systems, long propagation delays and rapid channel variation can render channel state information (CSI) outdated, thus requiring hybrid beamforming and resource allocation based on channel-aging [15]. Time-varying satellite visibility creates heterogeneous service conditions, with remote-region users experiencing fewer visible satellites, lower elevation angles, and longer service discontinuities than urban mid-latitude users [14]. By contrast, simultaneous multi-satellite visibility enables cooperative transmission and reception, motivating dynamic resource-allocation frameworks. These frameworks jointly optimize user scheduling, transmission power, combined weights, and bandwidth [16]. Traffic offloading is another established strategy for the management of NTN-TN resources, explored in the literature through joint interference management and coordinated coexistence [18] and traffic-aware scheduling with beam selection to increase throughput [19], collectively improving energy efficiency, capacity, coverage, and operational flexibility [17]. Interference Management for NTN–TN: While the integration of NTN nodes into terrestrial cellular systems can improve spectrum utilization and extend coverage beyond terrestrial infrastructure, it also creates complex interference conditions due to the distinct propagation characteristics of terrestrial and non-terrestrial links [29]. Several approaches have been proposed to manage such interference. In [20], a stochasticgeometry-based framework optimizes exclusion-zone radii for terrestrial cellular and backhaul networks coexisting with fixed-satellite services at 28 GHz. Alternatively, interference from terrestrial BS sidelobes can be mitigated by steering spatial nulls towards visible satellites [21]. More recently, [25] optimized downlink beamforming to protect incumbent VSATs while maximizing terrestrial sum rate, showing that aggressive spatial nulling can unnecessarily sacrifice capacity; however, the formulation does not consider interference coordination across geographically separated BSs. In [22], interference nulling and power control were considered for FR3 spectrum sharing, where terrestrial downlink transmissions interfere with LEO satellite uplinks. In [23], a joint detection, channel-estimation, and interference-cancellation framework was developed for terrestrial–satellite downlink coexistence. Although these CSI-dependent approaches can reduce interference at NTN receivers, they require accurate interference-

3

TABLE I S UMMARY OF PREVIOUS STUDIES ON RESOURCE AND INTERFERENCE MANAGEMENT IN INTEGRATED NTN–TN SYSTEMS . Reference [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] DT-STABLE

Interference Management ✗ ✗ ✗ ✗ ✓ ✗ ✓ ✓ ✓ ✓ ✗ ✓ ✓

FR2 band ✗ ✗ ✗ ✗ ✗ ✗ ✓ ✗ ✗ ✗ ✗ ✗ ✓

Environment Awareness ✓ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✓ ✗ ✓

Resource Variables

Objective

Satellite visibility, roaming CSI prediction Scheduling, power, bandwidth Cell switching, offloading Offloading, power, RBG Bandwidth, association, power Exclusion-zone radii, frequency allocation Beamforming, spatial nulling Power control, beamforming Beamforming, channel estimation Beam hopping, power allocation Beamforming Beam selection, association, scheduling

Access inequality analysis Mitigate channel aging Maximize sum rate Minimize power consumption Maximize spectral efficiency Maximize SLT, minimize energy Maximize terrestrial data rates Mitigate uplink interference Maximize power under INR/rate limits Maximize TN quality, limit interference Maximize throughput, fairness Maximize TN sum rate under INR constraints Maximize capacity under INR constraints

link CSI, increasing acquisition and processing overhead as the number of NTN receivers grows. Digital Twin for NTNs: The dynamic nature of NTN environments, characterized by time-varying satellite visibility, beam coverage, and interference conditions, requires accurate real-time network representation. High-fidelity simulators, such as the one previously proposed in [30], enable largescale evaluation of LEO constellations but remain largely limited to offline analysis. DTN bridge this gap through realtime bidirectional interaction between physical and virtual networks, thus enabling network mapping and closed-loop decision making [31], [32]. Recent advances include DTassisted collaborative resource allocation for multi-satellite LEO systems through coordinated beam hopping and power allocation [24], as well as DT-supported satellite-terrestrial edge computing for improved resource sharing, reduced computation delay, and lower satellite energy consumption [33]. Uniqueness of the Proposed Work: Table I compares the present research with previous studies, most of which addressed only a subset of the resource-management aspects relevant to integrated NTN–TN systems. While some prior studies incorporate site-specific environmental information, its role is largely confined to channel modeling and performance evaluation rather than resource-control decisions. Motivated by these research gaps, we develop an O-RAN-enabled, DT-assisted interference-management framework for contextaware spectrum sharing in the FR2 band. Unlike prior studies, the proposed framework integrates multi-domain contextawareness, proactive interference management at satellite receivers through DT-based prediction and what-if analysis, and closed-loop resource coordination over standardized ORAN interfaces. Collectively, these capabilities set apart the proposed framework from those previously proposed in the literature. III. S YSTEM M ODEL A. System Overview We consider an NTN–TN coexistence scenario, where terrestrial downlink transmissions share frequency resources with an uplink LEO satellite system. The terrestrial segment follows the O-RAN architecture, comprising multiple open radio units (O-RUs) with antenna arrays, connected to open distributed unit (O-DU) and open central unit (O-CU) entities

[34]. Radio resource control is supported by the near-realtime RIC (near-RT RIC), where xApps perform time-sensitive optimization, whereas the non-real-time RAN intelligent controller (non-RT RIC), hosted within the service management and orchestration (SMO) framework, provides non-real-time analytics, policy support, and DT-assisted contextual information. Let J = {1, 2, . . . , J} denote the set of O-RUs deployed over the considered geographic area. Furthermore, let U = {1, 2, . . . , U } denote the set of single-antenna user equipments (UEs). Each O-RU j ∈ J is equipped with a uniform rectangular array (URA) and a predefined discrete Fourier transform (DFT) beam codebook Bj = {1, 2, . . . , K} of candidate precoding vectors. Due to finite RF-chain and baseband processing resources, each O-RU can simultaneously activate only a limited number of beams [35]. The satellite segment consists of a set M = {1, 2, . . . , M } of visible LEO satellites operating over the same spectrum. Each satellite m ∈ M is modeled as a victim receiver that must be protected from harmful aggregate interference generated by terrestrial downlink transmissions. Therefore, the coexistence scenario requires coordinated O-RU beam activation and UE association decisions that enhance terrestrial service quality while protecting the coexisting satellite segment from harmful interference. In the proposed architecture, the CDT is maintained through the non-RT RIC using environmental, network, and satellitecontext information. It provides two context-aware CKMs for near-RT control: a 3D geo-environment aware terrestrial coverage map storing the received power associated with each beam over predefined spatial cells, and an interference map that quantifies the coupling of each O-RU–beam pair with visible satellite receivers. During near-RT operation, each UE position is mapped to the corresponding coverage-map cell, and the stored received power is used to derive the achievable rate for candidate O-RU–beam pairs. Concurrently, the perceived interference from each O-RU-beam pair to visible satellites is obtained from the precomputed interference map. A purely measurement-driven approach would require O(U JK) and O(JKM ) measurements per scheduling epoch. The proposed CDT instead reuses the precomputed maps and triggers only localized corrections when prediction mismatches occur, facilitating near-RT scalable operation. Assumptions: We adopt the following practical assump-

4

Fig. 1: O-RAN-enabled DT-STABLE framework for multi-timescale context-aware NTN–TN coexistence.

tions. A1: UE position is available to the DT-assisted O-RAN control plane. This is consistent with emerging 6G positioning capabilities, where UE locations can be obtained through multi-modal sensing [36] or GNSS and network-assisted fusion [37]. This is so because future 6G systems are expected to jointly support communication, localization, and sensing. This information is used to query DT-generated coverage maps and support context-aware BS–beam selection. Sensitivity of DT-STABLE to localization errors is evaluated in Section VI. A2: The DT can accurately estimate the interference generated by each BS beam towards LEO satellite receivers. Residual errors may arise from propagation-model mismatch, environmental changes, or antenna-pattern uncertainty. This uncertainty is assumed to follow a lognormal distribution in the linear domain, equivalently represented by a zero-mean Gaussian perturbation in the dB domain. Accordingly, the dB dB actual interference is modeled as Iactual = IDT + ϵI , where 2 ϵI ∼ N (0, σI ) and where σI denotes the standard deviation of the interference-prediction error in dB.1 The impact of DT interference-prediction errors is also evaluated in Section VI. A3: Inter-cell interference is neglected, which is reasonable for suburban macro deployments of downtilted BSs with typical inter-site distances and orthogonal frequency allocation among neighbouring cells via standard frequency reuse pattern [18]. Intra-cell inter-beam interference among simultaneously active beams is also neglected, as motivated by the structural orthogonality of the DFT codebook. However, in practice, residual inter-beam interference can arise due to multipathinduced leakage, which will be considered in our future work. B. DT-STABLE Framework Fig. 1 shows the architecture of the O-RAN-enabled, DTassisted framework for interference-aware NTN–TN coexistence. Within the non-RT RIC, the radio-twin rApp constructs and maintains long-term DT maps using slowly varying contextual information. The DT-generated CKMs are periodically refined by identifying prediction mismatch with RAN measurements collected by the near-RT RIC. The DTSTABLE xApp then combines the updated CKMs with the near-RT RAN state to jointly optimize beam activation and UE association under interference constraints. This separation 1 The Gaussian DT-error model is inspired by [38, Eq. (12)].

enables long-term map construction, mid-term feedback-driven refinement, and short-term resource optimization. A Walk-Through of DT-STABLE: The operation of the proposed framework is shown in Fig. 1. First, network and environmental information are collected through the O1 interface and stored in the context database at the SMO/nonRT RIC layer (Step 1). This information includes deployment maps, BS locations, beam codebooks, satellite ephemerides, weather information, and interference-management policies. Next, the radio twin rApp retrieves this information through the R1 interface and generates the channel knowledge map C and the interference map G (Step 2). The generated DT maps are then delivered to the near-RT RIC through the A1 interface (Step 3). At the near-RT RIC, the DT-STABLE xApp uses these maps together with RAN feedback received over the E2 interface, including UE positions, key performance measurements (KPMs), and measured rates. When a mismatch between DT-predicted and measured performance is detected, the xApp refines the DT information using limited real-world measurements before scheduling. Based on the updated DT and RAN information, the DT-STABLE xApp determines the beam activation and UE association decisions (Step 4). These decisions are sent to the relevant E2 nodes through the E2 interface and applied to the physical RAN (Step 5). The flow of Steps 1–5 establishes a closed-loop and proactive operation that synchronizes the digital and physical networks across multiple timescales. IV. DT-BASED C HANNEL AND I NTERFERENCE M ODELING IN NTN-TN C OEXISTENCE N ETWORKS A. Modeling of the Context-Aware Digital Twin (CDT) DT offers a practical approach for reducing dependence on extensive empirical measurements. This is done by synthesizing site-specific radio information that captures the underlying physical propagation environment [39], [40]. In the present study, the CDT is modeled using the following components. i) Real-world data: The CDT relies on contextual information from the real-world to represent the site. In this article, this information includes site geometry, terrain and building data, BS/O-RU locations, UE locations, satellite visibility, antenna configurations, and beam codebooks. These inputs allow the CDT to characterize both terrestrial service links and interference paths towards the satellite uplink. ii) Physics:

5

The physics layer models the site-specific TN and NTN propagation conditions relevant to coexistence control. Terrestrial O-RU coverage is generated using Sionna RT [27] over a high-fidelity virtual 3D environment,2 spatial coverage cells, antenna configurations, and beamforming patterns. The resulting multipath components are combined to obtain channel responses associated with the considered beams over the coverage region. During near-RT operation, UE locations are used only to identify the corresponding coverage-map cells. For the NTN links, O-RU–satellite RF signal propagation is modeled using OpenNTN [28]. Satellite positions are obtained from TLE data, while link losses are determined from the antenna responses, elevation angle, and atmospheric attenuation. iii) Network simulator: Instead of emulating the full network stack, the proposed CDT adopts a goal-specific link-level abstraction for NTN–TN coexistence control. The terrestrial propagation outputs are used to obtain the received-power information over the spatial coverage cells, while the NTN propagation outputs determine interference coupling between each O-RU beam and visible satellite receiver. These quantities are stored in terrestrial coverage map C and satellite interference map G, respectively, and are then used to compute the performance metrics, such as downlink throughput and aggregated interference. iv) Physical-to-virtual feedback: The physicalto-virtual feedback part transfers selected RAN observations back to the CDT. In the proposed O-RAN architecture, the near-RT RIC obtains UE locations, KPM reports, and measured link rates through the E2 interface. These observations are then compared to the CDT-predicted rates computed from the current channel map C. When the predictions deviate from the physical RAN observations, relevant link entries are refined using the latest feedback. This creates a midtimescale synchronization loop between physical RAN and the CDT, avoiding full map regeneration at every scheduling interval. v) DT-assisted resource control: The DT-STABLE xApp uses the current CDT state together with near-RT RAN feedback to determine beam activation and UE association decisions under the satellite interference constraints. The resulting control actions are then delivered to the physical RAN through the O-RAN control loop. Each terrestrial BS is equipped with an Nx × Ny URA and uses a finite DFT codebook for directional transmission. The antenna elements are separated by half-wavelength spacing, i.e., dx = dy = λ/2. Each antenna element follows the 3GPP TR 38.901 element pattern [8]. For a direction specified by elevation angle θ and azimuth angle ϕ, the element gain is modeled as shown below: ( 12



ϕ ϕ3dB

2

 + 12

θ − 90◦ θ3dB

Accordingly, the directional gain of the k-th terrestrial beam of BS j is computed as follows: 3GPP j,k GTN j,k (θ, ϕ) = Gelem (θ, ϕ) + Garray (θ, ϕ).

)

2 , Am

(1)

where Gmax is the maximum element gain, ϕ3dB and θ3dB are the horizontal while vertical half-power beamwidths, and Am is the maximum attenuation. 2 In this study, the virtual environment is constructed in Blender, incorporating BS/O-RU locations from OpenCell3D, building geometry from OpenStreetMap, and material properties assigned according to ITU standards.

(3)

For the BS–satellite direction (θj,m (t), ϕj,m (t)), the corresponding beam-dependent gain is computed as shown below:
TN GTN (4) j,k,m (t) = Gj,k θj,m (t), ϕj,m (t) . C. NTN Modelling for Coexistence Analysis 1) NTN Antenna Modelling: The receive antenna of each LEO satellite is modeled using a circular-aperture reflector pattern [1]. Let ψj,m (t) denote the off-axis angle between the boresight of the receive antenna of satellite m and the direction of the signal arriving from BS j. The corresponding normalized antenna power pattern is given by Eq. (5):  1,      J (k a sin(ψ 2 1 0 m j,m (t))) Fm (ψj,m (t)) = 4 ,  k0 am sin(ψj,m (t))     0,

ψj,m (t) = 0, 0 < |ψj,m (t)| ≤ 90◦ , otherwise, (5)

where J1 (·) is the first-order Bessel function of the first kind, am is the antenna-aperture radius of satellite m, while k0 = 2πfc /c is the wave number, with fc and c denoting the carrier frequency and the speed of light, respectively. The pattern is normalized such that Fm (0) = 1 along the antenna boresight. Assuming that the receiver noise figure and antenna noise temperature are independent of the signal arrival direction, the directional G/T of satellite m towards BS j is modeled as follows: 

B. TN Antenna Modelling for Coexistence Analysis

G3GPP elem (θ, ϕ) = Gmax − min

The k-th entry of BS j is represented by a unit-norm DFT precoding vector fj,k ∈ CNx Ny ×1 . Let aj (θ, ϕ) ∈ CM L×1 denote the URA steering response in direction (θ, ϕ), including the antenna element locations and BS orientation. The array gain induced by the k-th precoding vector is given by the following:   2 . (2) aH Gj,k j (θ, ϕ)fj,k array (θ, ϕ) = 10 log10

G T



 (t) = j,m

G T

 + 10 log10 Fm (ψj,m (t)) ,

(6)

m,max

where (G/T )m,max denotes the maximum boresight gain-tonoise-temperature ratio of satellite m. 2) NTN Channel Model: The NTN propagation channel is generated offline using OpenNTN [28], an open-source implementation of the 3GPP TR 38.811 geometry-based stochastic channel model (GSCM) [1]. For each BS–satellite link and ephemeris snapshot, channel generation is conditioned on the corresponding link geometry and deployment scenario. The GSCM determines the LOS/NLOS propagation state and generates the associated large-scale channel parameters, including shadow fading, delay and angular spreads, and the Rician K-factor. Conditioned on these parameters, the smallscale channel realization is synthesized from multipath clusters comprising rays with prescribed delays, powers, and angles of arrival and departure. The temporal evolution associated with

6

Fig. 2: Comparison of CKM construction approaches: (a) environment-aware multi-beam CKM; (b) environment-unaware CKM (with omnidirectional beam); and (c) environment-unaware multi-beam CKM.

satellite mobility is incorporated through the time-varying link geometry and the resulting Doppler characteristics. For each BS–satellite pair (j, m), the link geometry is determined from the BS position qj ∈ R3 and the timevarying satellite position qm (t) ∈ R3 , expressed in a common three-dimensional coordinate system. The satellite position is obtained from its ephemeris data, and the corresponding BS– satellite slant distance is dj,m (t) = ∥qm (t) − qj ∥2 . The BS–satellite geometry is further characterized by the elevation angle with respect to the local vertical at the BS. Let nj denote the local vertical unit vector at BS j. The elevation angle of satellite m observed from BS j is given by ! T (qm (t) − qj ) nj −1 . (7) θj,m (t) = sin dj,m (t) For directional antenna modeling, the BS–satellite direction is expressed in the local coordinate frame of the BS antenna panel, yielding the zenith and azimuth angles ϑloc j,m (t) and ϕloc (t), respectively. These angles account for the panel j,m orientation and downtilt. Let fj,k ∈ CNant ×1 denote the unitnorm precoding vector of beam k at BS j, and let aj (ϑ, ϕ) denote the corresponding array response, normalized such that ∥aj (ϑ, ϕ)∥22 = Nant . The directional transmit gain of beam k toward satellite m is then given by
loc loc Gtx j,k,m (t) = GE,j ϑj,m (t), ϕj,m (t) (8)
2 loc loc + 10 log10 aH j ϑj,m (t), ϕj,m (t) fj,k , where GE,j (ϑ, ϕ) denotes the antenna-element gain in dBi. Thus, (8) accounts for both the element radiation pattern and the beam-dependent array factor in the satellite direction. The large-scale propagation loss between BS j and satellite m is modeled as

the scenario- and elevation-dependent LOS probability. Conditioned on χj,m (t), the shadow-fading term is modeled as
2 SFj,m (t) ∼ N 0, σSF (E, θj,m (t), χj,m (t)) , (11) where E denotes the deployment scenario. The corresponding clutter loss is ( 0, χj,m (t) = LOS, CLj,m (t) = (12) CLE (θj,m (t), fc ) , χj,m (t) = NLOS, where CLE (·) denotes the scenario-dependent NLOS clutter loss. The atmospheric component P Lgj,m (t) accounts for gaseous, rain, and cloud attenuation along the Earth–space path according to the adopted 3GPP and ITU-R propagation models [1], [41]. For the considered configuration, the applicable scintillation term is set according to the corresponding 3GPP/ITU-R model. For each BS–satellite geometry, the GSCM further generates NMC independent stochastic channel realizations, with the link geometry held fixed. These realizations account for the smallscale multipath structure through the corresponding cluster and ray parameters, including delays, powers, angular characteristics, cross-polarization ratios, and random phases. Let (r) ℓj,k,m (t) denote the effective power coupling between beam k of BS j and satellite m obtained from the rth realization. The NTN coupling retained for interference evaluation is defined as follows: NMC 1 X (r) ℓ (t), NMC = 50. (13) ℓNTN j,k,m (t) = NMC r=1 j,k,m D. Construction of Multi-Antenna CKM

Let A denote the considered coverage region, which is parc titioned into Nc non-overlapping rectangular cells {Qn }N n=1 of size ∆x × ∆y , where ∆x = ∆y = 5 m in the considered implementation. Let Nc = {1, . . . , Nc } denote the g b s P LNTN (9) corresponding cell-index set. Let Bj = {1, . . . , Kj } denote the j,m (t) = P Lj,m (t) + P Lj,m (t) + P Lj,m (t), Nant ×1 is where P Lbj,m (t) denotes the basic path loss, P Lgj,m (t) denotes beam-index set of O-RU/gNB site j, where fj,k ∈ C the precoding vector associated with beam k ∈ B . Each Os j atmospheric attenuation, and P Lj,m (t) denotes the applicable scintillation loss. Building-entry loss is omitted since the RU/gNB site is equipped with an Nx × Ny URA, comprising Nant = Nx Ny antenna elements, and uses a finite DFT terrestrial BSs are deployed outdoors. beamforming codebook. The k-th precoding vector at site j is The basic path loss is expressed as √ expressed as fj,k = (ax (p) ⊗ ay (q))/ Nant , where ax (p) and P Lbj,m (t) = 32.45 + 20 log10 (fc [MHz]) + 20 log10 (dj,m (t)[km])ay (q) denote the horizontal and vertical DFT steering vectors, respectively; furthermore, (p, q) is the corresponding two+ SFj,m (t) + CLj,m (t), (10) dimensional DFT index pair, while ⊗ denotes the Kronecker where the first three terms correspond to the free-space path product. For the full DFT codebook, Kj = Nx Ny = Nant . loss. Let χj,m (t) ∈ {LOS, NLOS} denote the propagation We perform extensive site-specific high-fidelity RT to gencondition of the BS–satellite link, determined according to erate an offline coverage map for every candidate beam at each

7

O-RU/gNB site. For cell n, site j, and beam k, the average beamformed path gain is given Z by the following: 1 2 gn,j,k = hH (s)fj,k ds, (14) |Qn | Qn j where |Qn | denotes the area of cell Qn , s = (x, y) ∈ Qn denotes a spatial location within the cell, ds = dx dy is the infinitesimal surface element, and hj (s) denotes the channel vector between site j and, finally, location s. The resulting CKM enables DT-STABLE to proactively evaluate terrestrial service-side utilities, such as the received signal-to-noise ratio (SNR) and achievable spectral efficiency, for each candidate O-RU–beam–UE association. In this study, the complete multi-antenna CKM is represented by the beam-dependent received signal strength (RSS) tensor C = {ρn,j,k | n ∈ Nc , j ∈ J , k ∈ Bj }. where ρn,j,k = Pjb gn,j,k denotes the RSS in the linear domain for beam k of site j over cell n, while Pjb is the fixed transmit power allocated to one beam. During near-RT operation, the UE position is mapped to the corresponding coverage-map cell, and the stored RSS associated with each beam is used to approximate the received power at the UE location. The accuracy of this approximation depends on the cell size and the spatial variation of the received power within the cell. Fig. 2 compares three CKM models with varying levels of environmental and beam awareness. The environment-aware multi-beam CKM captures site-specific propagation conditions for each spatial cell and candidate beam. By contrast, the environment-unaware CKM relies only on location, and its multi-beam variant includes beam dependence without real propagation effects, thus generating idealistic uniform beam coverage. By incorporating site-specific information (blockage, terrain, and multi-path propagation), environment-aware multi-beam CKM can provide spatially consistent, physicsaccurate channel gain at any UE position within the O-RU coverage for each beam. E. Interference Map Construction We model the aggregate interference generated by active terrestrial BS beams toward each visible LEO satellite receiver. Let Pjb denote the fixed transmit power allocated to each active beam of BS j. Combining the effective NTN power coupling in Eq. (13) with the directional satellite gain-to-noisetemperature ratio in Eq. (6), the INR contribution of beam k at BS j toward satellite m is expressed in decibels as   dB Ij,k,m (t) = Pjb,dBW + 10 log10 ℓNTN j,k,m (t)
+ G/T j,m (t) − 10 log10 (κB).

(15)

where
Pjb,dBW is the per-beam transmit power in dBW, G/T j,m (t) is expressed in dB/K, B is the bandwidth in Hz, and κ is the Boltzmann constant. Accordingly, linear INR coupling coefficient is defined as follows: dB

gj,k,m (t) = 10Ij,k,m (t)/10 =

Pjb ℓNTN j,k,m (t) 10(G/T )j,m (t)/10 , (16) κB

where Pjb is expressed in watts and ℓNTN j,k,m (t) is the effective NTN power coupling defined in Eq. (13). Thus, gj,k,m (t) represents the linear INR contribution of beam k of BS j toward satellite m.

Cluster c = C, delay τC scatterer

LOS

O-RU (Tx)

LEO (Rx) RC rays

R1 rays

path

Cluster c = 1, delay τ1

(a) Geometry & Large-Scale Parameters

Small-Scale Cluster/Ray Parameters

Interference Map Construction

Cross-Polarization & Random Phases

Multipath Coeffs.

Array Response & Precoder

Random Coupling of Rays

Effective NTN Coupling

Directional Tx Gain

Departure / Arrival Angles

Per-Beam INR

LoS / NLoS Condition

Cluster Powers

Linear INR Coupling

Large-Scale Losses

Cluster Delays

Interference Map

BS j / Sat. m Geometry

FSPL, SF, clutter, atm.

(b) Fig. 3: NTN interference-map construction: (a) BS-to-LEO uplink channel model; (b) geometry, channel generation, and interference-map pipeline.

Next, let sj,k (t) ∈ {0, 1} indicate whether beam k of BS j is active. Assuming that all active beams overlap the considered bandwidth, transmit with unit activity, and contribute noncoherently at the satellite receiver, the aggregate INR at satellite m is obtained by summing the individual contributions in the linear domain as follows: X X Im (t) = sj,k (t)gj,k,m (t). (17) j∈J k∈Bj

The interference map provides a predictive representation of the interference coupling between terrestrial BS beams and visible satellite receivers over the satellite trajectory. By combining TLE-based satellite geometry with the directional antenna responses and NTN propagation model, the map captures the time-varying interference exposure associated with each candidate BS–beam pair. Fig. 3 illustrates the corresponding BS-to-LEO channel representation and the overall geometry-, channel-, and interference-map construction pipeline. This information enables DT-STABLE to anticipate interference events that may violate the satellite protection threshold and account for satellite protection in the near-RT resource-allocation policy. V. R ADIO R ESOURCE M ANAGEMENT FOR H ARMONIOUS NTN–TN C OEXISTENCE A. Optimization Problem Formulation In this section, we formulate the joint beam activation and UE association problem underpinning the proposed NTN–TN coexistence framework. The formulation balances terrestrial service performance and satellite protection, while promoting efficient use of available beam resources.

8

Decision Variables: We define binary beam activation variable sj,b ∈ {0, 1}, for all j ∈ J and b ∈ Bj , where sj,b = 1 indicates that beam b of BS j is active. We also define binary association variable xu,j,b ∈ {0, 1}, for all u ∈ U, j ∈ J , and b ∈ Bj , where xu,j,b = 1 indicates that UE u is served by BS j through beam P b. The load of beam b at BS j is defined as Kj,b = u∈U xu,j,b , for all j ∈ J and b ∈ Bj . Next, we define the achievable downlink rate for each candidate UE–BS–beam association. Let Pjb denote the transmit power allocated to one active beam of BS j. Let nu denote the coverage-map cell associated with UE u. The received power for candidate BS–beam pair (j, k) is obtained from the CKM as ρnu ,j,k . Accordingly, the received SNR is given by γu,j,k = ρnu ,j,k /σu2 , where σu2 = N0 BFN is the receiver noise power, N0 is the thermal noise power spectral density, B is the system bandwidth, and, finally, FN is the receiver noise factor. The corresponding achievable spectral efficiency is cu,j,k = log2 (1 + γu,j,k ) bps/Hz. To capture system throughput and user-level fairness, we adopt a logarithmic utility formulation [42]. We first define the loadaware terrestrial network utility as follows:   XX X cu,j,b . (18) U (x) ≜ xu,j,b log Kj,b u∈U j∈J b∈Bj

To capture the trade-off between terrestrial network utility and satellite interference protection, the joint optimization problem is formulated as shown in Eq. (19): X   th P0: max U (x) − λmax,m Im − Im + x,s

s.t.

m∈M X X   C1: xu,j,b = 1, ∀u ∈ U,    j∈J  b∈B j     C2: xu,j,b ≤ sj,b , ∀u, j, b,    X  C3: sj,b ≤ NRF , ∀j ∈ J ,  b∈Bj   X     C4: K = xu,j,b , ∀j, b, j,b    u∈U    C5: xu,j,b , sj,b ∈ {0, 1}, ∀u, j, b.

Algorithm 1: DT-STABLE Optimization Input: Candidate beam sets {Fj }j∈J ; DT-predicted rates {cu,j,b }; DT-predicted interference couplings {gj,b,m }; interference thresholds th {Im }; interference-price caps {λmax,m }; RF-chain limit NRF ; step-size parameters α0 , β, c, γ, αmax ; tolerance ε; maximum number of iterations Tmax . ⋆ Output: Beam activation {s⋆ j,b } and user association {xu,j,b }. (0) (0) 1 Initialization: Set λm = 0, ∀m ∈ M, and µj,b = 0, ∀j ∈ J , b ∈ Fj ; set t = 0 and ∆ = ∞. 2 while ∆ > ε and t < Tmax do (t) 3 Set xu,j,b = 0, ∀u ∈ U , j ∈ J , b ∈ Fj . 4 foreach u ∈ U do 5 Select the tentative BS–beam pair n o (t) ⋆ ⋆ (ju , bu ) ∈ arg max log(cu,j,b ) − µj,b . j∈J b∈Fj

(t) Set xu,j ⋆ ,b⋆ = 1. u u

6 7 8 9 10

end (t) Set sj,b = 0, ∀j ∈ J , b ∈ Fj . foreach j ∈ J do foreach b ∈ Fj do Compute the price-adjusted beam score i X (t) h X (t) (t) (t) xu,j,b log(cu,j,b ) − µj,b − λm gj,b,m . Qj,b = u∈U

11 12

(t)

Qj,b values among those satisfying X (t) xu,j,b > 0. u∈U 13 14 15

16 17 18

(t)

(t)

Set sj,b = 1, ∀b ∈ Aj . end Define the globally active BS–beam set n o (t) (t) A = (j, b) : j ∈ J , b ∈ Aj . foreach u ∈ U assigned to some (j, b) ∈ / A(t) do (t) Set xu,j,b = 0 for the discarded assignment. Select the globally active BS–beam pair n o (t) ′ ′ (ju , bu ) ∈ arg max log(cu,j,b ) − µj,b . (j,b)∈A(t)

(t)

19

Set xu,j ′ ,b′ = 1.

20

end foreach m ∈ M do Compute the aggregate interference and residual X X (t) (t) (t) (t) th Im = sj,b gj,b,m , vm = I m − I m .

u

21 22

(19)

u

j∈J b∈Fj 23

24

25 26 27 28

where [z]+ ≜ max(z, 0) and Im denotes the aggregate linear-domain INR at satellite m, as defined in (17). The coefficient λmax,m ≥ 0 controls the trade-off between terrestrial network utility and interference protection for satellite m. th Therefore, no interference penalty is incurred when Im ≤ Im , whereas violations above the threshold are increasingly penalized as λmax,m increases. Constraint C1 ensures that each UE is served by exactly one BS–beam pair. Furthermore, Constraint C2 restricts user association to active beams, i.e., xu,j,b = 1 only if sj,b = 1. Next, Constraint C3 limits the number of simultaneously active beams in each BS to NRF , the number of available RF chains per BS. In its turn, Constraint C4 defines the load Kj,b of each pair of BS–beams, as used in the load-aware utility. Finally, Constraint C5 enforces the binary nature of the beam activation and user association variables. Theorem 1: Problem P0 is NP-hard.

m∈M

end (t) Select Aj as the set of at most NRF beams in Fj having the largest

Compute the adaptive step size      (t) α0 (t) γ max{0,vm } αm = min 1+c e −1 , αmax . t+1 Update the interference price n n oo (t+1) (t) (t) (t) λm = min λmax,m , max 0, λm + αm vm . end foreach j ∈ J do foreach b ∈ Fj do Compute the beam demand and price-induced supply n o X (t) dem,(t) sup,(t) (t) Kj,b = xu,j,b , Kj,b = min |U |, exp(µj,b − 1) . u∈U

29

Update the load price (t+1)

µj,b 30 31 32

(t)

= µj,b −

i β h sup,(t) dem,(t) Kj,b − Kj,b . t+1

end end Compute the price-variation convergence metric ∆= λ

(t+1)

(t)

(t+1)

−λ

+ µ 2

33 34 35 36

−µ

(t)

. 2

Set t ← t + 1. end (t−1) (t−1) Set s⋆ and x⋆ j,b = sj,b u,j,b = xu,j,b . ⋆ return {s⋆ }, {x }. j,b u,j,b

Proof: Introduce an auxiliary variable tm ≥ 0 for each th satellite such that Im ≤ Im +tm . For λmax,m > 0, maximizing the objective forces tm to its smallest feasible value, i.e., th tm = [Im − Im ]+ . Hence, P the objective of P0 can be equivalently written as U (x) − m∈M λmax,m tm . We reduce the 0–1 knapsack problem to a restricted instance of P0.

9

Consider n items with values vj > 0, integer weights aj ≥ 0, and capacity W . Construct P0 with |J | = |U| = n, |M| = 1, (1) NRF = 1, and two beams per BS j: an item beam bj (0) and a dummy beam bj . Set cj,j,b(1) = evj , cj,j,b(0) = 1, j

j

and ci,j,b = δ = e−nV −1 for i ̸= j, where V = maxj vj . Moreover, set gj,b(1) ,1 = aj , gj,b(0) ,1 = 0, I1th = W , and j j λmax,1 = nV + 1. Assigning each UE to its corresponding dummy beam is feasible and yields objective zero. Any cross-BS association contributes at most log δ = −nV − 1, while the maximum total positive utility is at most nV ; hence, any solution containing a cross-BS association has negative objective and cannot be optimal. Therefore, each UE is associated with its corresponding BS, and the selected beam has unit load, i.e., Kj,b = 1. Let S denote the set of P BSs activating their item beams. The objective then becomes j∈S vj − λmax,1 t1 , P subject to j∈S aj ≤ W + t1 . Since the item weights aj and W are integers, any capacity violation implies t1 ≥ 1 and therefore yields an objective at most nV − (nV + 1) < 0, whereas the all-dummy solution achieves zero. Thus,Pt1 = 0 at the optimum, and the problem reduces to maxS j∈S vj P subject to j∈S aj ≤ W , which is exactly the 0–1 knapsack problem. Thus, since the construction is polynomial and the 0–1 knapsack problem is NP-hard, P0 is NP-hard. ■ B. DT-STABLE Algorithm Development Since P0 is computationally intractable due to NP-hardness, we propose a dual-level solution to P0 to obtain a suboptimal yet efficient solution. The detailed derivation of the twolevel solution is provided in Appendix A. Based on this, Algorithm 1 proceeds as follows: • Lines 1–2: Initialization. Interference prices λm and load prices µj,b are initialized to zero, together with the iteration index and convergence metric. • Lines 3–6: Provisional UE association. Given the current load prices, each UE selects the candidate BS–beam pair with the highest load-adjusted rate. This provisional association provides an estimate of the demand and terrestrial utility associated with each candidate beam before enforcing the RF-chain constraint. • Lines 7–13: Beam scoring and activation. Each BS evaluates its candidate beams using the load-adjusted utilities of their provisionally associated UEs and the corresponding satellite-interference price. The occupied beams are ranked according to their resulting scores and, up to NRF , beams with the largest scores are activated. • Lines 14–23: UE reassignment. UEs provisionally associated with deactivated beams are reassigned to the globally active BS–beam pair providing the highest load-adjusted rate. Since the candidate beams in this step are already active, the reassignment does not change the aggregate satellite interference. • Lines 24–26: Interference-price update. For each satellite, the aggregate interference generated by active beams is compared with the corresponding interference threshold. A positive interference residual increases the corresponding price, thereby imposing a larger penalty on highly interfering

beams in the next beam-selection step, whereas a negative residual decreases the price. The adaptive update strength is controlled by α0 , c, γ, and αmax , while the projection bounds the interference price within its prescribed interval. • Lines 27–31: Load-price update. For each BS–beam pair, the associated UE demand is compared with its priceinduced supply. When the demand exceeds the supply, the load price increases, reducing attractiveness of the congested beam in the next association step. Otherwise, the load price decreases to reflect the lower congestion level. • Lines 32–35: Convergence and output. The convergence metric is evaluated from the changes in interference and load prices. Iterations terminate when this metric falls below ε or when Tmax is reached. The resulting beam-activation and UE-association decisions are then returned. Computational Complexity of Algorithm 1: In each iteration, provisional UE association and beam-load evaluation require O(U JK) operations, while the interference-aware beam-score computation requires O(JKM ). Beam selection, interferenceprice updating, and load-price updating require O(JK log K), O(JKM ), and O(U JK) operations, respectively. Hence, over at most Tmax iterations, the overall complexity is O(Tmax (U JK + JKM + JK log K)). Since P0 is a nonconvex mixed-integer problem, global optimality of Algorithm 1 is not guaranteed. Instead, DTSTABLE provides a computationally efficient heuristic that adaptively balances terrestrial capacity and satellite interference protection, particularly when protecting additional satellite links from harmful interference would require a substantial sacrifice in network utility. SectionVI-G evaluates the runtime and optimization performance of DT-STABLE against the Gurobi benchmark. C. O-RAN Compliant Multi-Time Scale Optimization The proposed DT-assisted control framework operates across multiple O-RAN timescales to support context-aware NTN–TN coexistence. The non-RT RIC maintains the longterm DT model of the radio environment and uses slowly varying contextual information to construct and update terrestrial CKMs and satellite interference maps. The near-RT RIC, by contrast, processes faster network dynamics, including UE locations, serving-cell information, and measured KPMs collected from the RAN through the E2 interface [13]. This separation enables the DT to combine long-term environmental information with near-RT feedback from the physical network. The DT maintains two CKMs (C and G) to predict terrestrial beam coverage and interference at satellite receivers. These maps are updated by the non-RT RIC using site-specific context and satellite information, and the resulting feasible beam sets are transferred to the near-RT RIC through the A1 interface. The near-RT xApp then combines the DT predictions with E2 KPM feedback to perform closed-loop resource control. A key objective of the proposed multi-timescale design is to avoid exhaustive beam sweeping. Conventional beam sweeping may require testing a large number of BS–beam pairs, which increases measurement overhead and control latency [45]. By contrast, the DT provides a ranked prediction

10

Algorithm 2: Multi-Timescale O-RAN DT Resource Optimization Input: Beam codebooks {Bj }j∈J DT coverage map C DT interference map G LEO TLE data TLEO context state η interference thresholds th {Im } top-Kp probing size Kp mismatch threshold δmis measured KPM feedback RF-chain limit NRF . ⋆ Output: Optimized beam activation {s⋆ j,b } and user association {xu,j,b }. 1 Initialization: Set Fj ← Bj , ∀j ∈ J . Umis ← ∅ and Ku ← ∅, ∀u ∈ U . 2 Long-term DT map generation at the Non-RT RIC: 3 Collect context information through O1 and update η. 4 Update satellite visibility using TLEO and obtain the visible satellite set M. 5 Update the DT coverage map C = {ρn,j,b }n∈Nc , j∈J , b∈F based on the j current context state η. 6 Update the DT interference map G = {gj,b,m }j∈J , b∈F , m∈M based on j satellite visibility and BS beam patterns. 7 Provide C, G, and {Fj } from the Non-RT RIC to the Near-RT RIC via the A1 interface. 8 Mid-term DT correction and selective probing at the Near-RT RIC: 9 Collect UE locations, serving-cell information, and measured rate KPMs {r̂u }u∈U via the E2 interface. 10 foreach u ∈ U do 11 Map UE u to its corresponding coverage-map cell nu . 12 Initialize the UE-specific DT received-power predictions as ρ̃u,j,b ← ρnu ,j,b , ∀j ∈ J and ∀b ∈ Fj . 13 Compute the corresponding DT-predicted rates   ρ̃u,j,b , ∀j ∈ J , ∀b ∈ Fj . cu,j,b ← log2 1 + 2 σu cur cur 14 Obtain the current serving pair (ju , bu ) and the corresponding cur ,bcur . DT-predicted rate r̃u ← cu,ju u 15 Compute the relative prediction error |r̂u − r̃u | , eu ← max(r̃u , ϵr ) where ϵr > 0 is a small numerical constant. 16 if eu > δmis then 17 Umis ← Umis ∪ {u}. 18 Select the top-Kp candidate BS–beam pairs according to the DT-predicted received power Ku ← TopKj∈J , b∈Fj (ρ̃u,j,b , Kp ) . 19 20 21 22

23

24 25 26 27 28

end end Probe the candidate sets {Ku }u∈Umis through E2 and obtain the measured received powers {ρ̄u,j,b }(j,b)∈Ku . Replace the corresponding DT-predicted received-power estimates for the current scheduling instance, ρ̃u,j,b ← ρ̄u,j,b , ∀u ∈ Umis and ∀(j, b) ∈ Ku . Recompute the corresponding achievable rates   ρ̃u,j,b , ∀u ∈ Umis , ∀(j, b) ∈ Ku . cu,j,b ← log2 1 + 2 σu Short-term scheduling at the Near-RT RIC: th Run Algorithm 1 with inputs {Fj }, {cu,j,b }, {gj,b,m }, {Im }, U , M, and NRF . ⋆ ⋆ Obtain {sj,b } and {xu,j,b }. ⋆ Dispatch {s⋆ j,b } and {xu,j,b } from the Near-RT RIC to the BSs via the E2 interface. ⋆ return {s⋆ j,b }, {xu,j,b }.

of promising BS–beam pairs before measurement. When a UE’s measured rate is consistent with the DT prediction, the current DT map is reused without triggering additional probing. Beam probing is performed only for UEs whose measured KPMs exhibit a predefined deviation from the DTpredicted KPMs due to residual DT–real-world mismatch. For these UEs, the near-RT xApp selects only the top-Kp feasible BS–beam pairs predicted by the DT and refines the corresponding map entries using selective measurements. This DT-guided beam-sweeping procedure reduces the search space while maintaining adaptability to mobility, blockage, and environmental changes. The main steps of Algorithm 2 are summarized as follows: • Lines 1–4: Initialization and long-term DT update. The non-RT RIC collects context information through O1, updates the DT context state, and determines the visible satellite set from the LEO TLE data. • Lines 5–6: DT map update. The DT updates the coverage

TABLE II S IMULATION , CHANNEL - MODELING , AND DT-STABLE PARAMETERS . Parameter

Value Network Parameters Terrestrial BSs / UEs 19 ( [43]) / 50–300 LEO satellites / elevation range 20 / 25.01◦ –32.46◦ Deployment area 4 × 4 km2 Carrier frequency / bandwidth 28 GHz / 50 MHz [3] Per-beam transmit power / UE noise figure 33 dBm / 7 dB Monte Carlo runs / snapshots / interval 150 / 10 / 1 s BS / UE height 25 m / 1.5 m UE mobility / speed Random waypoint / 0.1–2 m/s Terrestrial BS Antenna Parameters BS array / element spacing 4 × 4 URA / 0.5λ (H/V) Beam codebook / RF chains 16 2-D DFT beams / 4 Element pattern / downtilt 3GPP TR 38.901 [8] / 12◦ Element gain / 3-dB beamwidth 8 dBi / 65◦ (H/V) UE antenna Single-element dipole Sionna RT and CKM Parameters Ray-tracing platform Sionna RT [27] Ray samples / maximum path depth 5 × 106 / 5 Diffraction / CKM resolution Enabled / 5 × 5 m2 Material models ITU-R P.2040 brick, wood, and concrete NTN Channel and Interference Parameters Constellation / ephemeris Starlink / CelesTrak TLE [44] Minimum visibility elevation 25◦ NTN model / scenario Urban [1] Channel realizations per link 50 Clutter loss (LOS / NLOS) 0 dB / elevation- and frequencydependent Atmospheric attenuation Gaseous, rain, and cloud Additional losses Scintillation / building entry / misc.: 0 dB Satellite G/T / INR threshold 13 dB/K / −6 dB† DT-STABLE and Multi-Timescale Parameters Step sizes (α0 , β) (0.001, 0.001) Price parameters (γ, c) (5.0, 1.2) Caps (αmax , λmax,m ) (0.08, 50) Maximum iterations / tolerance 300 / 10−3 Minimum UE spectral efficiency 0.5 bit/s/Hz Mismatch threshold / probing size δmis = 0.3 / Kp = 5 † The −6 dB threshold corresponds to an INR of 25%, consistent with

commonly adopted in satellite coexistence studies [21].

map C = {ρn,j,b } using the current context state and updates the interference map G = {gj,b,m } based on the visible satellites and BS beam patterns. • Lines 7–9: A1 and E2 information exchange. The nonRT RIC provides the updated DT maps and beam sets to the near-RT RIC through A1. The near-RT RIC collects UE locations, serving-cell information, and measured rate KPMs through E2. • Lines 10–20: DT mismatch detection and candidate selection. For each UE, the xApp maps its location to the corresponding coverage-map cell and retrieves the DTpredicted received power for each candidate BS–beam pair. The achievable rates are then computed. The measured rate of the current serving pair is compared with its DT-predicted rate, and if the relative error exceeds δmis , the UE is added to Umis . Then, the top-Kp candidate BS–beam pairs are selected according to their DT-predicted received powers. • Lines 21–23: Selective probing and DT prediction correction. For UEs in Umis , only the selected candidate BS–beam pairs are probed. The measured received powers replace the corresponding DT-predicted received-power estimates for the current scheduling instance, and the associated achievable rates are recomputed. • Lines 24–28: DT-STABLE scheduling and control. The resulting achievable rates, interference couplings, feasible

11

TABLE III O PTIMIZATION TRADE - OFF AND NORMALIZED BENCHMARK PERFORMANCE . Approach Interference-Optimal DT-STABLE Rate-Optimal

INR Sat. (%) 99.99 89.63 66.99

Nominal BS Downtilted BS

7.5 5.0 2.5 0.0 25-45

(a) Elevation-angle distribution of visible satellites.

(a) Optimization trade-off statistics Median INR (dB) 95th INR (dB) -14.35 -7.94 -13.09 -4.56 -9.64 1.84

45-70

70-90

(b) Impact of MBS antenna sidelobes on the worst-satellite INR.

Fig. 4: Satellite visibility and sidelobe-induced interference characteristics in the considered coexistence scenario. 1.0 0.8

0.8

0.6

CDF

CDF

1.0

Int.-Opt. DT-STABLE Rate-Opt.

0.4 0.2 0.0

Int.-Opt. DT-STABLE Rate-Opt.

0.6 0.4 0.2

0

100

200

300

Total Sum Rate [bps/Hz]

(a)

Mean Rate (bps/Hz) 133.55 210.96 257.24

(b) Ideal-Point Score Comparison (U = 100) DT-STABLE NL-Mod. NL-Str. 0.740 0.695 0.615 0.951 0.970 1.000 0.813 0.784 0.728

5G Baseline 1.000 0.587 0.708

Worst-satellite INR p95 [dB]

Metric Norm. Rate Norm. INR Sat. Norm. Ideal-Point Score

Mean INR (dB) -15.14 -13.07 -9.16

0.0 −30

−20

−10

0

10

INR [dB]

(b)

Fig. 5: CDF comparison of the benchmark schemes: (a) terrestrial sum rate and (b) satellite interference.

beam sets, and INR thresholds are provided to the DTSTABLE xApp, which executes Algorithm 1. The resulting beam-activation and UE-association decisions are sent to the BSs through the E2 interface. VI. P ERFORMANCE E VALUATION A. Simulation Setting The proposed DT-STABLE framework was evaluated using a synchronized CDT–physical-network simulation environment. The CDT provided site-specific radio information to support decision making, while the physical-network layer evaluated the resulting beam activation and UE association decisions. Robustness of DT-STABLE against practical impairments is evaluated in Section VI-F. The simulation, channelmodeling, Sionna RT/CKM, and DT-STABLE multi-timescale control parameters are summarized in Table II. We considered NTN–TN spectrum coexistence, where terrestrial downlink transmissions share resources with a LEO satellite system. The site-specific environment was modeled as a DT of a 4×4 km2 region near Québec City, Canada. The TN comprised a total of 19 macro O-RUs/MBSs whose geographic coordinates were obtained from the OpenCellID database [43]. Following the 3GPP suburban macro assumptions, each macro O-RU had a height of 25 m and a 1000 m service region. Each O-RU was equipped with a half-wavelength-spaced directional 4×4 URA and a finite DFT codebook of 16 candidate beams.

PBS-L 0.703 0.739 0.720

95th Rate (bps/Hz) 206.40 273.47 304.24 PBS-M 0.672 0.867 0.749

PBS-H 0.638 0.984 0.744

For the considered LEO constellation, satellite positions were derived from TLE data obtained from CelesTrak [44]. Satellite visibility across the deployment area was determined using a minimum elevation threshold of 25◦ . Under this condition, Fig. 4a shows that the elevation-angle distribution of the visible Starlink satellites is concentrated toward the lower-elevation range. Fig. 4b further examines the impact of elevation angle on terrestrial interference by incorporating the MBS antenna response into the calculation of the worst-case INR. While antenna downtilt reduced the p95 INR by up to 5.19 dB at elevations above 70◦ , the reduction decreased to 3.09 dB in the 25◦ –45◦ range, leaving the residual p95 INR above the −6 dB protection threshold. Thus, low-elevation satellites constitute the most critical coexistence regime in the considered scenario, as they are both more prevalent in the visibility set and less effectively protected by antenna downtilt. To provide a conservative worst-case assessment of DT-STABLE, the ensuing analysis considered visible satellites with the lowest elevation angles at each observation. B. Optimization Strategies By tuning the adaptive interference-price parameters (γ, c, αmax ) in Algorithm 1, which control how strongly the interference-prices respond to constraint violations, the proposed algorithm realizes three operating points: interference optimal, balanced, and rate optimal. The interference-optimal point prioritizes satellite protection, whereas the rate-optimal point prioritizes terrestrial throughput. Fig. 5(b), Fig. 5(a), and Table III(a) characterize the rate–interference trade-off among the three operating points. Fig. 5(b) presents the INR CDF, where the vertical dashed line marks the target threshold of −6 dB, while Table III(a) summarizes the corresponding INR satisfaction and terrestrial sum-rate performance. For operating scheme a, theP INRPsatisfaction ratio was defined as ηaINR = PM S R (r,s) 100% m=1 1[INRa,m ≤ INRth ], where R, S, r=1 s=1 RSM and M denoted the numbers of simulation runs, snapshots, and evaluated satellites, respectively. The interference-optimal configuration achieved the highest INR satisfaction, 99.99%, but the lowest mean system sum-rate, 133.55 bps/Hz. In other words, nearly every evaluated satellite remained below the INR threshold under this configuration. By contrast, the rateoptimal configuration achieved the highest mean system sumrate, 257.24 bps/Hz, while its INR satisfaction decreased to 66.99%. DT-STABLE provided a balanced operating point, achieving 89.63% INR satisfaction and 210.96 bps/Hz. Compared to the rate-optimal setting, it improved INR satisfaction by 22.64 percentage points, while increasing the sum-rate

5 0 −5 −10 −15 −20 −25 −30

INR Satisfaction (%)

Sum Rate (bps/Hz)

INR (dB)

12

350 300 250 200 150 100

. . . L 5G T-ST od -Str BS- BS-M BS-H P P D NL-M NL P

(a)

30 25 20

50

100

150

200

250

Number of Users (U)

300

70 60 50

50

80 70 60 50 40 100

150

200

250

Number of Users (U)

300

Fig. 8: Average user-rate loss vs. the number of users. Fig. 9: QoS satisfaction vs. the number of users.

by approximately 1.58× relative to the interference-optimal setting. The balanced point, referred to as DT-STABLE, was used throughout the remainder of the evaluation, as it provides an effective balance between terrestrial throughput and satellite interference protection. C. Benchmark Schemes We compared DT-STABLE with the three benchmark schemes in terms of satellite INR, terrestrial sum rate, average user-rate loss, and QoS satisfaction. To capture different levels of INR-constraint enforcement, the no-load-aware scheme was evaluated under moderate and strict filtering levels, while the per-BS beam selection scheme was assessed under low, medium, and high local beam-penalty levels. • Baseline 1: 5G Baseline. This scheme follows the conventional 3GPP beam management procedure, where each BS exhaustively sweeps all candidate DFT beams and each user is assigned to the BS–beam pair with the highest channel gain. Satellite interference was not considered, making this a standard-compliant, but coexistence-unaware benchmark that served as the terrestrial-rate upper bound. • Baseline 2: No-Load-Aware (NL) Scheme. This benchmark does not account for load-dependent bandwidth sharing among users served by the same BS–beam pair. Thus, multiple users can be assigned to the same beam without an explicit load penalty. Two variants, moderate and strict, were considered by changing the weight of the satellite INR term in the optimizer. • Baseline 3: Per-BS Beam (PBS) Selection. This benchmark performs interference-aware beam selection locally at each BS, inspired by [21]. PBS-Low, PBS-Medium, and PBS-High denote low, medium, and high local beampenalty levels, respectively, with higher levels providing more conservative satellite protection. Unlike DT-STABLE, PBS does not coordinate aggregate interference across BSs or balance the load across active beams. Fig. 6(a) presents the aggregate INR distribution for U = 100. The 5G baseline yielded a mean INR of −6.53 dB with only 55.32% satisfaction, showing that rate-driven allocation

150

200

250

300

Number of Users (U) NL-Mod. NL-Str.

PBS-L PBS-M

PBS-H

Fig. 7: INR satisfaction vs. the number of users.

90

50

100 5G DT-ST.

Sum Rate (bps/Hz)

35

80

(b)

QoS Satisfaction (%)

Avg. Rate Loss (%)

40

90

. . . L 5G T-ST od -Str BS- BS-M BS-H P P D NL-M NL P

Fig. 6: Benchmark distributions: (a) satellite INR and (b) terrestrial sum rate. 45

100

400

DT-STABLE

3GPP

50

150

300 200 100 0 100

200

250

300

Number of Users (U)

Fig. 10: DT-STABLE vs. 3GPP sum rate.

alone cannot ensure NTN–TN coexistence. DT-STABLE reduced the mean INR to −13.07 dB with 89.63% satisfaction, achieving INR control comparable to the NL schemes while avoiding excessive beam suppression. By contrast, PBS-Low under-controls INR, whereas PBS-High provides more conservative interference protection. Fig. 6(b) shows the corresponding sum-rate distribution. Although the 5G baseline achieved the highest rate because it is unconstrained by INR, DT-STABLE reached 210.96 bps/Hz, outperforming NL-Moderate/Strict, PBS-Medium/High, and the other schemes providing comparable INR protection. This rate–interference trade-off is also reflected in Table III(b). For each scheme a, the rate and INR-satisfaction performance were normalized by their respective maxima as Ranorm = R̄a / maxq∈A R̄q and Panorm = ηaINR / maxq∈A ηqINR , so that both ratios are dimensionless and (1, 1) represents the theoretically perfect solution. We then score Sa = 1 − p defined the normalized ideal-point √ (1 − Ranorm )2 + (1 − Panorm )2 / 2, which measures each scheme’s Euclidean distance to this ideal point, with Sa ∈ [0, 1] and larger values indicating closer proximity. DTSTABLE attained the highest score, Sa = 0.813, indicating the closest proximity to the ideal point among the evaluated benchmarks. Fig. 7 shows a similar trend across user densities. The 5G baseline achieved approximately 55% INR satisfaction for all U , while the PBS schemes improved satisfaction through higher penalty levels at the expense of capacity. The NL schemes also became slightly less effective as the number of users increased; for example, the moderate setting decreased from 91.65% to 88.43%. DT-STABLE maintained a more balanced behavior, with INR satisfaction decreasing from 91.09% at U = 50 to 87.56% at U = 300, as it jointly incorporated DT-predicted rates, INR prices, and load prices. Fig. 8 further illustrates this trade-off. The more conservative NL schemes provided stronger interference protection but incurred larger user-rate loss relative to the 5G baseline at lower user densities, whereas the PBS schemes exhibited

Total Swept Beams

Sum rate (bits/s/Hz)

13

250 200 Static DT map DT-STABLE

150 0

20

Blockage window Trigger

40

60

Time (slot)

80

100

Fig. 11: Sum-rate comparison during temporary blockage.

increasing rate loss as the number of users increased. At U = 100, DT-STABLE achieved 89.63% INR satisfaction with a 25.99% rate loss relative to the 5G baseline and maintained higher INR satisfaction than PBS-Low across all evaluated user densities. Fig. 9 reports QoS satisfaction across user densities, defined as the percentage of UEs whose achieved spectral efficiency meets the minimum requirement rmin = 0.5 bit/s/Hz. DTSTABLE’s QoS satisfaction decreased from 73.91% to 49.73% as U increased from 50 to 300, while remaining more stable than the other INR-aware schemes. At U = 100, it achieved 68.67% QoS satisfaction, outperforming the PBS schemes and remaining competitive with the NL schemes. These results indicate that load pricing helps limit excessive user concentration on active BS–beam pairs. D. Advantages of Digital Twin Predicted Channel We evaluated the effect of the DT-predicted, site-specific channel model by comparing it to a conventional 3GPP urban macrocell (UMa) baseline. The DT channel maps were generated using Sionna RT in a reconstructed 3D environment, incorporating actual BS locations, terrain and building geometry, material-dependent propagation, and multi-antenna beam patterns. This approach enabled the DT model to capture channel variations due to blockage, reflection, diffraction, and scattering. Fig. 10 compares the total sum rate achieved with the DT-predicted channel and the statistical 3GPP baseline at different user densities. The DT-based model consistently outperformed the 3GPP baseline. With an increase of the number of users from 50 to 300, the average sum rate with the DTpredicted channel rose from 146.85 to 315.64 bps/Hz, whereas the 3GPP baseline increased only from 25.49 to 43.45 bps/Hz. The median sum rate followed a similar trend, with the DTbased model achieving 147.07–318.07 bps/Hz as compared to 24.47–42.87 bps/Hz for the 3GPP abstraction. Moreover, the relative gain of the DT-based model widened from approximately 5.8× to 7.3× as user density increased. These results demonstrated that channel fidelity is essential for sitespecific resource allocation. Unlike statistical models, the DTbased channel model captured physics-grounded, deploymentspecific RF conditions and provided link-quality estimates that closely reflected the actual BS–UE environment, thus enabling more effective beam selection and user association. E. Advantages of Multi-time Scale Digital Twin Assisted Resource Management To assess DT-assisted beam management under short-term radio variations, we consider a 120-slot time-domain experiment based on the 28-GHz DT environment. The channel state is held fixed over each 10-slot feedback interval and

TABLE IV R ECOVERED RATE VERSUS BEAM SWEEPING BUDGET.

105 104 103 102

Full Sweep

Selective DT-STABLE

Swept beams, Kp Normalized rate 1 0.670 2 0.793 3 0.841 5 0.867 8 0.919 16 1.000

Fig. 12: Total number of swept beams.

updated at the beginning of the next interval, consistent with the adopted DT-refinement timescale. During slots 60–90, a temporary mixed-blockage event is introduced for 30% of the UEs. Among the affected UEs, a subset experiences complete blockage of the scheduled links, whereas the remaining links undergo partial attenuation. This creates a controlled mismatch between the offline CKM and instantaneous radio conditions. The static DT continues to rely on the offline CKM, whereas DT-STABLE uses measurement feedback to identify and correct inaccurate link predictions. At each feedback instant, the DT-predicted rate of the serving link is compared with the UE measurement. When the relative prediction error exceeds the predefined mismatch threshold, DT-STABLE triggers selective probing of DT-ranked candidate beams and updates the corresponding CKM information before resource allocation. Fig. 11 shows the temporal sum-rate evolution of the static DT and DT-STABLE under the blockage experiment. The setup considers 50 UEs, 19 BSs, 16 candidate beams per BS, and feedback evaluation every 10 slots. Prior to the blockage interval, the two schemes exhibit comparable performance because the offline CKM remains consistent with the channel conditions. Once the blockage is introduced, the static DT cannot account for the resulting channel mismatch and continues to make resource-allocation decisions from outdated link information. In contrast, DT-STABLE incorporates the feedback-driven CKM refinement in Algorithm 2 and adapts the affected beam and association decisions accordingly. Over the blockage interval, the average sum rate increases from 172.40 bits/s/Hz with the static DT to 197.91 bits/s/Hz with DT-STABLE, corresponding to a 14.8% improvement. Fig. 12 quantifies the corresponding beam-measurement overhead. A full sweeping strategy evaluates all BS–beam candidates for every UE at each feedback instant, requiring 182,400 measurements. Selective exhaustive sweeping limits these measurements to links with a detected prediction mismatch, reducing the total to 32,528 measurements. DTSTABLE further restricts the search to the top-Kp candidates identified from the DT. With Kp = 5, the 107 detected mismatch events require only 535 beam measurements, corresponding to reductions of 99.71% and 98.36% relative to full and selective exhaustive sweeping, respectively. Table IV further evaluates the quality of the DT-guided beam ranking under different probing budgets. For each mismatch event, the highest rate among the top-Kp CKM-ranked beams of the serving BS is normalized by the highest rate from exhaustive probing of all 16 beams, and the resulting ratios are averaged over all mismatch events. With Kp = 8, the DT-guided candidate set achieves 91.9% of the full-probing rate on

100

50 0

5

Positioning Error (%)

10

Per-satellite INR (dB)

150

QoS Satisfaction (%)

Sum Rate (bps/Hz)

14

90 80 70 60 0

5

Positioning Error (%)

(a)

10

(b) Pre-opt. DT-STABLE

8 6 4 2 0

Low

Mid

High

Elevation Profile

−10 −20 −30 0

1

2

3

5

100

150

Uncertainty σ (dB) Fig. 14: Impact of interference-map uncertainty on INR.

0

500

−10 −20 −30 −40

INR threshold

0

Runtime (ms)

10

INR (dB)

INR Exceedance (dB)

Fig. 13: Impact of UE positioning error: (a) sum rate and (b) QoS satisfaction.

10

INR threshold

Low Mid High 25–45° 45–70° 70–90°

400 300 200 100 0

50

75

Number of Users

Fig. 15: Mean INR exceedance vs. elevation profiles. Fig. 16: INR distribution across elevation profiles. Fig. 17: DT-STABLE runtime versus user density.

average, demonstrating that measurement overhead can be substantially reduced while retaining most of the rate available from exhaustive beam probing. F. Robustness Evaluation We further evaluated the robustness of DT-STABLE against practical impairments. In this evaluation, the scheduler operated with mismatched DT information, while the resulting beam activation and user association decisions were validated over the real-world layer. We examined the following two main sources of DT mismatch: errors in UE positioning and uncertainty in the interference map. We also analyzed how the satellite elevation angle influences interference and coexistence performance by altering the BS–satellite geometry. 1) Impact of Imperfect UE Positioning Estimation: Fig. 13 shows the impact of UE positioning errors on DT-STABLE performance for 50 users. Since the CKM was queried based on UE locations, an error in the estimated position may cause the scheduler to use channel information that is not fully aligned with the actual UE location. This can lead to less effective BS–beam associations and reduce the achieved terrestrial performance. With accurate UE positioning, DTSTABLE achieved a mean total sum rate of 147.33 bps/Hz and a QoS satisfaction ratio of 88.4%. When the positioning error increased to 5%, these values decreased to 87.88 bps/Hz and 75.8%, respectively. At 10% positioning error, the mean total sum rate and QoS satisfaction ratio further decreased to 64.06 bps/Hz and 68.6%, respectively. Overall, the total sum rate decreased by 56.5%, whereas the QoS satisfaction ratio decreased by 19.8 percentage points. This indicates that positioning errors primarily degraded the achievable rate, while their effect on the ability to satisfy the minimum userrate requirement was comparatively smaller. These results highlight the importance of accurate spatial information for the performance of the DT-assisted scheduler. 2) Impact of Interference Prediction Uncertainty: Fig. 14 evaluates DT-STABLE under the relaxed Assumption A2 for

100 users, where σI denotes the standard deviation of the prediction error applied to the DT beam-to-satellite coupling coefficients. The scheduler operates using the nominal interference map, while the actual couplings are modeled as true DT 10ϵj,k,m /10 , with ϵj,k,m ∼ N (0, σI2 ). Thus, gj,k,m = gj,k,m σI = 0 dB represents an exact interference map, while larger values introduce increasing coupling uncertainty. The terrestrial channel rates are kept unchanged to isolate this effect. As σI increases, the INR distribution shifts toward higher values. At σI = 0 dB, the mean and median INR are −13.01 dB and −13.05 dB, respectively. At σI = 5 dB, they increase to −10.88 dB and −10.92 dB, while the INR satisfaction ratio decreases from 89.66% to 79.15%, demonstrating the impact of interference-map uncertainty on satellite protection. 3) Impact of Elevation Angle: We next evaluated how the satellite elevation angle affects DT-STABLE’s worst-case INR performance. This is relevant because BS–satellite interference depends on the propagation geometry and the MBS antenna response toward the satellite direction. Fig. 15 shows that, before optimization, the worst-satellite INR excedence is elevation dependent. The worst-case INR for all three elevation profiles exceeded the target INR threshold, with mean exceedances of 7.2 dB, 5.1 dB, and 0.8 dB for the low- (25◦ –45◦ ), mid- (45◦ – 70◦ ), and high-elevation (70◦ –90◦ ) cases, respectively. After DT-STABLE scheduling, these exceedances were significantly reduced to 1.8 dB, 1.5 dB, and 0.2 dB, corresponding to reductions of 75.0%, 70.6%, and 75.0%, respectively. Fig. 16 further shows the post-scheduling per-satellite INR distributions for the three elevation profiles. The median INR remained comparable across the low-, mid-, and high-elevation cases, with values of −16.47 dB, −16.17 dB, and −17.55 dB, respectively. In contrast, the upper tail of the INR distribution improved with increasing elevation, as the 95th-percentile INR decreased from −4.54 dB at low elevation to −5.11 dB at mid elevation and −8.52 dB at high elevation. Accordingly, the INR satisfaction ratio increased from 89.63% at low elevation to 92.47% and 98.50% at mid and high elevation,

15

TABLE V P ERFORMANCE COMPARISON BETWEEN G UROBI AND DT-STABLE. #

BS

U

1 2 3 4 5 6 7 8 9 10

10 10 12 12 14 14 16 16 19 19

100 150 200 220 220 250 240 270 200 260

Runtime Gurobi (s) DT-STABLE (s) 3.264 0.207 5.395 0.348 13.815 0.484 15.157 0.543 15.386 0.601 21.214 0.640 11.558 0.653 14.646 0.697 11.493 0.594 16.861 0.757

SR/BS (bit/s/Hz) Gurobi DT-STABLE 15.23 15.29 13.81 14.65 13.38 14.09 13.00 12.60 13.26 12.60 13.13 12.36 12.66 12.87 12.38 12.59 15.34 13.92 13.30 14.05

INR Satisfaction Gurobi DT-STABLE 100% 75% 100% 80% 100% 75% 100% 85% 100% 70% 100% 80% 100% 70% 100% 75% 100% 70% 100% 75%

respectively, indicating that DT-STABLE maintained effective satellite protection across all elevation profiles, with lower upper-tail interference at higher elevation profiles. G. Computational Complexity and Gurobi Benchmarking Since DT-STABLE was executed by the near-RT resourceallocation xApp, its runtime is critical for practical deployment. Fig. 17 reports computation time for different UE populations. As expected, the runtime increased with the number of UEs because a larger user set expands the BS–beam association space. However, it remained below 0.45 s in all cases, including the 150-UE scenario, indicating compatibility with the near-RT RIC control timescale of 10 ms–1 s. To evaluate its computational efficiency and solution quality, DT-STABLE was compared with Gurobi on reduced instances of Problem P0. Table V reports the runtime, sum rate per BS, and INR satisfaction for different BS and UE configurations. Gurobi solves the mixed-integer problem using branch-andbound, which searches the combinatorial beam-activation and UE-association space while enforcing the constraints of Problem P0. Therefore, it achieved 100% INR satisfaction, but required significantly higher runtime. Gurobi required 3.26– 21.21 s, whereas DT-STABLE needed only 0.207–0.757 s. DT-STABLE avoids this combinatorial search through iterative price-based beam selection and user association, achieving much lower runtime with comparable sum-rate performance, although its INR satisfaction was up to 30 percentage points lower than that of Gurobi. For example, DT-STABLE achieved 14.09 versus 13.38 bps/Hz/BS for 12 BSs and 200 UEs, and 12.60 versus 13.00 bps/Hz/BS for 12 BSs and 220 UEs. In six of the ten configurations, DT-STABLE attains a slightly higher sum rate than Gurobi since it does not strictly enforce the INR constraints for all satellites. These results show that DTSTABLE provides a practically implementable solution within the near-RT time constraint, with limited sum-rate and INRsatisfaction degradation compared with exact optimization. VII. C ONCLUSION In this paper, we introduced DT-STABLE, an O-RANenabled DT framework for NTN–TN spectrum coexistence. Using site-specific channel and satellite-interference maps, DT-STABLE was found to support proactive beam selection and user association without exhaustive real-time measurements. We formulated the resource allocation problem with practical constraints, showed it to be NP-hard, and developed a low-complexity Lagrangian-inspired price-based algorithm for near-real-time control. Integrated into a multi-timescale O-RAN control loop, DT-STABLE combines long-term DT mapping, mid-term corrections, and short-term scheduling. Site-specific simulations over the 28 GHz band showed that

DT-STABLE achieved an effective balance between terrestrial throughput and satellite protection. For 100 users, DTSTABLE achieves 89.63% INR satisfaction and a mean sum rate of 210.96 bps/Hz, while the framework managed to attain nearly complete INR satisfaction when satellite protection was prioritized. Under dynamic blockage, DT-STABLE increased throughput by 14.8% as compared to a static DT map while reducing measurement overhead by 99.71%. It also satisfied near-RT RIC time budget with 15.5x to 33.1x faster execution than the exact optimization. Taken together, these results demonstrate that combining site-specific DT knowledge with O-RAN control enables scalable, efficient, and interferenceaware resource optimization for NTN–TN spectrum sharing. A PPENDIX A D ERIVATION OF DT-STABLE A. Lower-Level Resource Allocation For a fixed interference-price vector λ(t) , the term P (t) th m λm Im is constant with respect to (x, s). Hence, the lower-level problem is X max U (x) − λ(t) (20) m Im (s). x,s

m∈M

P Using Kj,b = u xu,j,b , the utility becomes X X U (x) = xu,j,b log(cu,j,b ) − Kj,b log Kj,b . u,j,b

(21)

j,b

Relaxing xu,j,b , sj,b ∈ {0, 1} to [0, 1] yields a convex program, since −K Plog K is concave and C1–C3 are affine. Dualizing Kj,b = u xu,j,b with multiplier µj,b ∈ R gives X L= xu,j,b [log(cu,j,b ) − µj,b ] u,j,b

+

X

[µj,b Kj,b − Kj,b log Kj,b ] −

j,b

X

λ(t) m sj,b gj,b,m .

m,j,b

(22) sup The stationary condition gives Kj,b = min{|U|, exp(µj,b − P dem 1)}. Comparing it with Kj,b = u xu,j,b gives the load-price update  β  sup,(t) dem,(t) (t+1) (t) . (23) Kj,b − Kj,b µj,b = µj,b − t+1 For fixed load prices, the association-dependent term is linear in xu,· over the association simplex. Thus, an optimum can be selected at a simplex vertex, yielding h i (t) (t) (ju⋆ , b⋆u ) = arg max log(cu,j,b ) − µj,b , xu,j ⋆ ,b⋆ = 1. u

j,b

u

(24) The corresponding interference-adjusted beam score is X (t) X (t) (t) Qj,b = xu,j,b [log(cu,j,b ) − µj,b ] − λ(t) m gj,b,m . (25) u

m (t) Qj,b

Beams with ≤ 0, or with no associated UE, are deactivated. If more than NRF beams remain at a BS, only the NRF largest-score beams are retained. For the fixed tentative association, this maximizes the retained beam-score sum subject to the RF-chain limit. UEs on discarded beams (t) are reassigned using (24) over A(t) = {(j, b) : sj,b = 1}, restoring C1–C3. We emphasize that this sequential restoration is polynomial-time but does not nessesarily guarantee global optimality, consistent with the NP-hardness of P0.

16

B. Upper-Level Adaptive Interference Pricing After the lower-level allocation, the upper level computes (t) (t) th vm = Im − Im and    i (t) α0 h (t) 1 + c eγ[vm ]+ − 1 , αmax , (26) αm = min t+1 followed by h i (t) (t) (t) λ(t+1) = Π λ + α v (27) [0,λmax,m ] m m m m . A positive residual increases λm , penalizing subsequent interference toward satellite m, whereas a negative residual decreases it. The projection keeps λm ∈ [0, λmax,m ]. The two levels iterate until convergence or t = Tmax , yielding a computationally tractable but not certifiably globally optimal solution to P0. R EFERENCES [1] 3GPP, “Study on New Radio (NR) to Support Non-Terrestrial Networks (Release 15),” 3rd Generation Partnership Project (3GPP), Technical Report TR 38.811, 2020. [Online]. Available: https://portal.3gpp.org [2] Samsung Research, “6G Spectrum: Expanding the Frontier,” Samsung Research, White Paper, May 2022. [3] 3GPP, “NR; base station (BS) radio transmission and reception,” 3GPP TS 38.104, Version 18.5.0, Release 18, 2024. [4] Federal Communications Commission, “SpaceX Gen2 NGSO satellite system, order and authorization,” FCC 22-91, 37 FCC Rcd 14882, Dec. 2022. [5] A. Yastrebova et al., “Theoretical and Simulation-Based Analysis of Terrestrial Interference to LEO Satellite Uplinks,” in Proc. IEEE GLOBECOM, 2020, pp. 1–6. [6] T. Miller et al., “Expanding digital connectivity through satellite broadband in the 28 ghz band,” Plum Consulting, Tech. Rep., Oct. 2021. [Online]. Available: https://plumconsulting.co.uk/expandingdigital-connectivity-through-satellite-broadband-in-the-28-GHz-band/ [7] R. C. Hansen, “Array pattern control and synthesis,” 2004. [Online]. Available: https://api.semanticscholar.org/CorpusID:2804502 [8] 3GPP, “Study on Channel Model for Frequencies from 0.5 to 100 GHz,” 3rd Generation Partnership Project (3GPP), Technical Report TR 38.901, 2020, Version 16.1.0. [9] A. Molisch and D. A. Barnhart, “NewSpectrum: Interference from terrestrial cellular systems to satellites: Measurement, modeling, and effective mitigation strategies,” NSF Award Number 2433981. Directorate for Engineering, Division of Electrical, Communications and Cyber Systems. 2025., p. 33981, Apr. 2025. [10] A. M. Seid et al., “A multi-agent DRL-based dynamic resource allocation in O-RAN-enabled TN-NTN metaverse services,” IEEE Trans. Commun., vol. 73, no. 12, pp. 14 243–14 259, 2025. [11] K. Alam et al., “A comprehensive tutorial and survey of O-RAN: Exploring slicing-aware architecture, deployment options, use cases, and challenges,” IEEE Commun. Surveys Tuts., vol. 28, pp. 1637–1678, 2026. [12] X. Lin et al., “6g digital twin networks: From theory to practice,” IEEE Communications Magazine, vol. 61, no. 11, pp. 72–78, 2023. [13] S. Seth et al., “Post: Programmable open spectrum twin for enhanced spectrum sharing in 6g bands,” IEEE Commun. Standards Mag., vol. 10, no. 3, pp. 37–45, 2026. [14] M. Almekhlafi, A. Lesage-Landry, and G. K. Kurt, “Access inequality in LEO satellite networks: A case study of high-latitude coverage in northern québec,” IEEE Open J. Veh. Technol., vol. 6, pp. 1613–1630, 2025. [15] O. Abbasi and G. Kaddoum, “Channel aging-aware LSTM-based channel prediction for satellite communications,” IEEE Netw. Lett., vol. 6, no. 3, pp. 183–187, 2024. [16] Q. Peng et al., “Dynamic resource allocation in distributed mimo-leo satellite networks,” arXiv preprint arXiv:2505.20891, 2025. [17] M. Ozturk, M. Salamatmoghadasi, and H. Yanikomeroglu, “Integrating terrestrial and non-terrestrial networks for sustainable 6G operations: A latency-aware multitier cell-switching approach,” IEEE Netw., vol. 40, no. 2, pp. 156–164, 2026. [18] M. M. Rahman et al., “Joint interference management and traffic offloading in integrated terrestrial and non-terrestrial networks,” IEEE Trans. Commun., vol. 73, no. 8, pp. 6748–6763, 2025.

[19] H. Alam et al., “Optimizing integrated terrestrial and non-terrestrial networks performance with traffic-aware resource management,” arXiv preprint arXiv:2410.06700, 2024. [20] A. U. Rahman, M. A. Kishk, and M.-S. Alouini, “Coexistence of terrestrial and satellite networks in the 28-ghz band,” IEEE Trans. Aerosp. Electron. Syst., vol. 59, no. 6, pp. 8342–8354, 2023. [21] S. Kang et al., “Terrestrial-satellite spectrum sharing in the upper midband with interference nulling,” in Proc. IEEE ICC. IEEE, 2024, pp. 5057–5062. [22] M. Tsampazi et al., “Satellite-terrestrial spectrum sharing in FR3 through QoS-aware power control and spatial nulling,” in Proc. IEEE Int. Symp. Spectrum Innovation (DySPAN), 2026, pp. 1–10. [23] S. Jia et al., “Joint detection, channel estimation and interference nulling for terrestrial-satellite downlink co-existence in the upper mid-band,” in Proc. IEEE GLOBECOM, 2025, pp. 2771–2776. [24] R. Zhao et al., “Demand-aware beam hopping and power allocation for load balancing in digital twin empowered LEO satellite networks,” IEEE Trans. Wireless Commun., vol. 24, no. 6, pp. 5084–5098, 2025. [25] K. Allouch et al., “Beamforming design for cellular coexistence with incumbent satellite terminals in the upper mid-band,” IEEE Open J. Commun. Soc., 2026. [26] Y. Zeng et al., “A tutorial on environment-aware communications via channel knowledge map for 6g,” IEEE Communications Surveys & Tutorials, vol. 26, no. 3, pp. 1478–1519, 2024. [27] J. Hoydis et al., “Sionna RT: Differentiable Ray Tracing for Radio Propagation Modeling,” in Proceedings of the IEEE Conference, Dec. 2023, pp. 317–321. [28] T. Düe et al., “OpenNTN: An open-source framework for non-terrestrial network channel simulations,” in Proc. 28th Int. Workshop Smart Antennas (WSA), 2025, pp. 1–7. [29] B. Shang et al., “Spectrum sharing in satellite-terrestrial integrated networks: Frameworks, approaches, and opportunities,” IEEE Netw., vol. 40, no. 4, pp. 162–171, 2026. [30] H. Zhou et al., “A high-fidelity and high-efficiency simulator for 6Gintegrated space–ground networks,” Engineering, vol. 56, pp. 62–78, 2026. [31] Y. Wu, K. Zhang, and Y. Zhang, “Digital twin networks: A survey,” IEEE Internet Things J., vol. 8, no. 18, pp. 13 789–13 804, 2021. [32] H. Gao et al., “Digital twin enabled 6g radio testing: Concepts, challenges and solutions,” IEEE communications Magazine, vol. 61, no. 11, pp. 88–94, 2023. [33] Z. Ji, S. Wu, and C. Jiang, “Cooperative multi-agent deep reinforcement learning for computation offloading in digital twin satellite edge networks,” IEEE J. Sel. Areas Commun., vol. 41, no. 11, pp. 3414–3429, 2023. [34] O-RAN Working Group 1, “O-RAN Architecture Description 5.00,” ORAN Alliance e.V., TS 5.00, 2021. [35] A. F. Molisch et al., “Hybrid beamforming for massive mimo: A survey,” IEEE Communications magazine, vol. 55, no. 9, pp. 134–141, 2017. [36] A. Alkhateeb et al., “Deepsense 6g: A large-scale real-world multimodal sensing and communication dataset,” IEEE Communications Magazine, vol. 61, no. 9, pp. 122–128, 2023. [37] M. Brambilla et al., “Integration of 5g and gnss technologies for enhanced positioning: An experimental study,” IEEE Open J. Commun. Soc., vol. 5, pp. 7197–7215, 2024. [38] A. Bazzi and M. Chafii, “How much training is needed with a digital twin?” arXiv preprint arXiv:2609.01220, 2026. [39] A. Alkhateeb, S. Jiang, and G. Charan, “Real-time digital twins: Vision and research directions for 6g and beyond,” IEEE communications magazine, vol. 61, no. 11, pp. 128–134, 2023. [40] H. Luo et al., “Digital twin aided massive MIMO CSI feedback: Exploring the impact of twinning fidelity,” IEEE Trans. Commun., vol. 73, no. 12, pp. 15 909–15 925, 2025. [41] International Telecommunication Union, “Propagation data and prediction methods required for the design of Earth–Space telecommunication systems,” International Telecommunication Union, Geneva, Switzerland, Recommendation ITU-R P.618-14, Aug. 2023. [42] Q. Ye et al., “User association for load balancing in heterogeneous cellular networks,” IEEE Trans. Wireless Commun., vol. 12, no. 6, pp. 2706–2716, 2013. [43] OpenCellID, “OpenCellID,” https://www.opencellid.org, 2025, accessed: 2025-07. [44] CelesTrak, “CelesTrak,” https://celestrak.org/, 2026, accessed: 2026-0528. [45] M. Qurratulain Khan et al., “Machine learning for millimeter wave and terahertz beam management: A survey and open challenges,” IEEE Access, vol. 11, pp. 11 880–11 902, 2023.

Related documents

• Multi-Timescale Context-Aware Radio Resource Optimization for Harmonious NTN-TN Coexistence: A Digital Twin and Open RAN Enabled Framework#1448343
• Rate-Splitting for Service-Aware Multi-Orbit Orchestration in 6G TN-NTN Networks#1084568
• Transport Layer Performance Combining Multipath and Network Coding in Mobile Satellite Networks#1215312
• LLM-Assisted Hierarchical Interference Control for Satellite-Terrestrial Integrated Networks#1297460
• Enhancing 5G NTN VSAT RACH in GNSS-Denied environments#1077093