F (l) = σ(W (l)F (l−1) + b(l)) (35) where F (l) is the output of the layer l, W (l) is the convolution kernel, b(l) is the bias term, and σ is the ReLU activation function. Temporal dependencies across frames are captured using 1D temporal convolution, as described by Eq. (36).
Tj = T ∑ t=1 Kj(t) · Xt (36) Where Kj(t) is the 1D kernel applied across the time axis, and Tj is the temporal feature vector. The reconstructed output ˆ X ∈RH×W representing the turbulence-compensated beam, is obtained using a transposed convolution (deconvolution) operation, as defined in Eq. (37).
ˆ X = Deconv(F (N)) (37) Where F (N) is the final feature map from the last convolutional block. The training objective minimizes the pixel-wise MSE between the predicted output and the clean target image Xclean, given by Eq. (38).
LMSE = 1 HW H ∑ i=1 W ∑ j=1 ( ˆ Xi,j −Xclean,i,j )2 (38) To preserve fine structural details, we include a spatial gradient loss, formulated in Eq. (39).
Lgrad = ∇ˆ X −∇Xclean 2 2 (39) The total loss function is then defined by Eq. (40)
Ltotal = λ1LMSE + λ2Lgrad (40) where λ1 and λ2 are hyperparameters determined empirically during training, optimization is performed using the Adam optimizer with a learning rate µ. The training convergence is monitored using the Root Mean Square Error (RMSE), calculated as shown in Eq. (41).
RMSE = 1 N N k=1 ˆYk −Yk 2 (41) In this equation ˆYk and Yk denote predicted and actual beam intensities for the training sample. Additionally, the Pearson correlation coefficient (CC) is used to assess the structural similarity, as expressed in Eq. (42)
CC = ∑ ( ˆ X −¯ X)(Xclean −¯ Xclean) √∑ ( ˆ X −¯ X)2 ∑ (Xclean −¯ Xclean)2 (42) (b) Time-domain correlation sequence generation model for turbulence fading channel To generate the time-domain correlation sequence model for a turbulence fading channel, we consider the large-scale Cln A(ρ) and small-scale Cln B(ρ) logarithmic covariance functions. The ACF model for optical turbulence can be expressed by Eq. (43)
Cintensity(ρ) = exp[Cln A(ρ) + Cln B(ρ)] −1 (43) where the Cintensity(ρ) represents the intensity of ACF at spatial separation ρ. The comprehensive expressions for Cln A(ρ) and Cln B(ρ) are provided in Eq. (44) and Eq. (46) through Eq. (48) and Eq. (49).
Cln A(ρ) = 1.06σ2 variance ˆ 1 0 ˆ ∞ 0 ζ−11/6 exp[−ζ ζA Z0(ρ√ωζ D )(1 −cos ζϕ)]dζdϕ (44)
Cln A(ρ) ≃0.16σ2 varianceζ7/6 A F1(7 6; 1; −ωρ2ζA 4D ) (45)
Cln B(ρ) = 1.06σ2 variance ˆ 1 0 ˆ ∞ 0 ζ−11/6 exp[−ζ ζB Z0(ρ√ωζ D )(1 −cos ζϕ)]dζdϕ (46) Scientific Reports | (2026) 16:8921 12 | https://doi.org/10.1038/s41598-026-40704-2 www.nature.com/scientificreports/
Cln B(ρ) ≃1.27σ2 variance ( ωρ2 DζB )5/12 K5/6( √ ωρ2ζB D ) (47) Where σ2 variance is the Rytov variance. D represents the atmospheric structure constant. ω = 2π/λ′ is the wave number, with λ′ being the wavelength. F1(.) signifies a confluent hypergeometric function. K5/6 is denoted the modified Bessel function of the second kind. The parameters ζA and ζB are given by:
ζA = 2.61 1 + 1.11σ12/5 variance (48)
ζB = 3 ( 1 + 0.69σ12/5 variance ) (49) By substituting Eq. (44)-(49) into Eq. (43) and applying the transformation given in Eq. (50),
ρ = W⊥τ ′ (50) where, W⊥ is the transverse wind speed and τ ′ is the time index, we can derive the time-dependent ACF as shown in Eq. (51).
Cintensity(τ ′) = exp [ 0.49σ2 variance (1 + 1.11σ12/5 variance)7/6 F1 ( 7 6; 1; −ωW 2 ⊥(τ ′)2ζA 4D ) + 0.50σ2 variance (1 + 0.69σ12/5 variance)5/6 ( ωW 2 ⊥(τ ′)2ζB 4D )5/12 K5/6 (√ ωW 2 ⊥(τ ′)2ζB 4D )] −1 (51) The normalized ACF Ψintensity(τ ′) is defined by Eq. (52).
Ψintensity(τ ′) = Cintensity(τ ′) Cintensity(0) (52) The normalized ACF represents a stationary stochastic process and its power spectrum (S) form can be derived as shown in Eq. (53).
Sintensity(Ω′) = 4 ˆ ∞ 0 Ψintensity(τ ′) cos(Ω′τ ′)dτ ′ (53) Here, Ω′ = 2πv′, where v′ = 1/λ′ denotes the light frequency. To further discretize the data, a normalized digital frequency ϕ = 2π/fsample is obtained by sampling at fsample. When added with Gaussian noise AWGN, at a certain power. σ2 noise passes through a rational filter transfer function H(ϕ) as given by Eq. (54). The temporal signal can then be generated using the relations in Eq. (55) and (56).
Sintensity(ϕ) = ∞ n=1 Ψintensity(n)e−jϕn = B(ϕ) A(ϕ) σ2 noise = H(ϕ)σ2 noise (54)
A(ϕ) = 1 + a1e−jϕ + · · · + ane−jnϕ (55)
B(ϕ) = 1 + b1e−jϕ + · · · + bme−jmϕ (56) From the covariance structure of the autoregressive moving average (ARMA) process, it can be deduced from Eq. (57).
Ψfilter(0) + N∑ n=1 afilterΨfilter(−n) = σ2 noise, Ψfilter(m) + N∑ n=1 anΨfilter(m −n) = 0 (57) Let the relations in Eqs. (58) and (59) hold.
αn′ = Ψintensity(n + 1) + Ψ∗ intensity,1×nϕn×1 (58)
[ Ψ(n+1)×(n+1) intensity Ψintensity(−n −1) Ψn×1 intensity Ψintensity(0) ] [ 1 ϕn×1 ]
[σ2 noise 0 ] (59) Using a recursive solution method, we can obtain the results shown in Eqs. (60), (61), and (62). Scientific Reports | (2026) 16:8921 13 | https://doi.org/10.1038/s41598-026-40704-2 www.nature.com/scientificreports/
ϕ(n+1)×1 = [ϕn 0 ]
- ξn+1 [ϕn 0 ] (60) σ2 noise(n + 1) = σ2 noise(n) ( 1 −| ξ(n + 1)|2) (61) σ2 noise(n)ϕ(n)= Z(n) (62) Where, ξ(n + 1) = −αn′ /σ2 noise(n) allows us to determine the coefficients, ϕn of the turbulence fading channels transfer function, and calculate the time-domain correlation sequence Z using the above given equations. However, this calculation is based on AWGN and does not account for the probability density function (PDF), meaning the generated sequence may not conform to turbulence disturbance theory. To address this, assuming a free-space receiver without aperture smoothing, the Gamma-Gamma (G-G) function is used as the PDF for the turbulence fading channel. Its parameters βGG and γGG can be related to the turbulence structure constant C2 n as shown in Eq. (63)-Eq. (66).
Pintensity (J)= 2(βGGγGG)(βGG +γGG )/2 Γ(βGG)Γ(γGG) J ( J ⟨J⟩ )(βGG +γGG )/2 (63)
KβGG−γGG (2 √ βGGγGGJ ⟨J⟩ ), J> 0 (64)
βGG = 1 σ2 A = 1 exp(σ2 ln A) − 1 (65)
γGG = 1 σ2 B = 1 exp(σ2 ln B ) − 1 (66) Where the expressions for the logarithmic variances σ2 ln A and σ2 ln B , which are required for these calculations, are provided in Eqs. (67) and (68).
σ2 ln A = 0.49σ2 variance ( 1+1 .11σ12/5 variance )7/6 (67)
σ2 ln B = 0.51σ2 variance ( 1+0 .69σ12/5 variance )9/6 (68) Here, ⟨Jn⟩ represents the nth moment of the light intensity J, and Γ(·) denotes the Gamma function. When the normalized random light intensity sequence { ˆJn/⟨ ˆJ⟩ } 1×N From the simulation, the G-G PDF follows, which the vector relation can express in Eq. (69).
ˆJ = ˆJ1/⟨ ˆJ⟩ ˆJ2/⟨ ˆJ⟩ ... ˆJn/⟨ ˆJ⟩ ∈ PD FGamma−−Gamma (69) Where i =1 , 2,...,n . By rearranging ˆJ according to the rank of Z, the CC of the sequence ρJ,Z can be calculated using the formula given in Eq (70).
CCρJ,Z = ⟨(J(i) −⟨ J⟩)(Z(i) −⟨ Z⟩)⟩√ ⟨(J(i) −⟨ J⟩)2(Z(i) −⟨ Z⟩)2⟩ =1 (70) This perfect mapping ensures that the temporal correlation information from Z is transferred to J, thereby satisfying the dual criteria defined in Eq. (71). ˆJ ∈ P DFGamma−−Gamma, Ψintensity (n)= ⟨(J(i) −⟨ J⟩)(Z(i − n) −⟨ Z⟩)⟩ (71) Eq. (51) verifies that the ACF generated by ˆJ, Ψintensity (n) matches the numerical simulation of the theoretical solution Eq. (52), while preserving the PDF information. Filter coefficients are obtained using the Y ule–W alker algorithm, and a temporal Gaussian correlation sequence is produced by integrating A WGN. Finally, through the rank mapping of the time-domain Gaussian correlation, the random light intensity sequence following the G–G PDF is reordered to accurately represent turbulence disturbance information across frequency bands. Scientific Reports | (2026) 16:8921 14| https://doi.org/10.1038/s41598-026-40704-2 www.nature.com/scientificreports/ Multiplexing and demultiplexing techniques To maximize the throughput of inter-satellite FSO links, we propose a hybrid approach that combines WDM with OAM-based MDM. This system is modeled to operate natively in the mid-infrared (mid-IR) spectral region, leveraging its advantages for FSO communications. The implemented model relies on compact, single- layer dielectric metasurface structures, specifically designed for the mid-IR, which are modeled to function as integrated transceivers for efficient OAM multiplexing and demultiplexing, generating and receiving orthogonal coaxial OAM beams as depicted in Fig. 3. (a) WDM in Mid-IR WDM enhances channel capacity by simultaneously transmitting multiple data streams on different optical wavelengths. For the mid-IR regime, this is efficiently achieved using native optical components which provide direct emission within the 3-5µm atmospheric transmission window. This approach avoids the inefficiencies of wavelength conversion and is well-suited for high-capacity, multi-gigabit-per-second data links. (b) OAM-based MDM in Mid-IR OAM-based MDM boosts capacity by using orthogonal spatial modes to transmit multiple channels with minimal crosstalk. In the mid-IR, OAM beams of different orders remain orthogonal. Spiral phase plates (SPPs) and compact near-IR OAM generation methods can be adapted for mid-IR use. Combining WDM with OAM- based MDM enables higher capacities through simultaneous transmission of multiple high-capacity channels across wavelengths and orthogonal OAM modes. (c) Optical Metasurface-based OAM Multiplexing/ Demultiplexing in Mid-IR To achieve a highly integrated and compact form factor for our FSO transceiver, we model single-layer dielectric metasurfaces designed to operate in the mid-infrared spectrum. These sub-wavelength structures provide unprecedented control over the phase, amplitude, and polarization of light, enabling them to efficiently generate and separate multiple coaxial OAM beams. For multiplexing, the metasurface is modeled to impart a spatially varying phase profile that converts an array of incident Gaussian beams (each potentially at a different WDM wavelength) into a set of coaxial OAM beams with distinct topological charges (l). The phase function for generating an OAM beam with charge l is given by ϕ(x, y) = l · arctan (y/x ), which defines the required spatial phase distribution. For demultiplexing at the receiver, the complementary metasurface model performs the inverse operation. It applies a phase correction that transforms each incoming OAM mode ℓ into a focused Gaussian spot at a distinct lateral position on a multi-pixel detector array. This spatial separation allows for the simultaneous detection of all multiplexed channels with minimal crosstalk. The metasurface’s efficiency in our model is characterized by its ability to direct optical power into the desired diffraction order. The performance is quantified by the Signal-to- Crosstalk Ratio (SXR) between different OAM modes, which we optimize to be >15 dB across the operational mid-IR band (e.g., 3-5 µm). Unlike bulky conventional optics like spiral phase plates, the modeled metasurfaces represent an ultra-thin (on the order of the wavelength) and lightweight alternative, suggesting a potential path for future compact, next-generation FSO terminals that are suitable for size, weight, and power (SWaP)-constrained applications like satellite communications. To contextualize the trade-offs between conventional SLM-based approaches and the proposed metasurface implementation, Table 3 compares their key performance parameters. This comparison highlights the advantages of metasurfaces for compact, SWaP-constrained applications despite their static phase profiles. Both approaches present practical implementation considerations that would require careful engineering in physical prototypes, but are suitable for the system-level modeling presented in this work. Fig. 3. Multiplexing and Demultiplexing.
Scientific Reports | (2026) 16:8921 15 | https://doi.org/10.1038/s41598-026-40704-2 www.nature.com/scientificreports/