Wavefront Astigmatism and Coherent-Aperture Limits in HF Skywave Propagation: An IONORT Case Study

Abstract

This paper reports a simulation study of the two-dimensional (2D) phase structure of high-frequency (HF) skywave wavefronts arriving at ground level after single-hop oblique ionospheric reflection, bounding the aperture over which a single-radius spherical wavefront model remains consistent for distributed long-baseline ground-based HF receiver arrays. For a representative short-range link (transmitter at Beijing; ground range 419 km; 10 MHz; IRI-2020 epoch 2026-04-14T06:00 UTC), three-dimensional IONORT ray tracing produces a dense landing-phase footprint fitted to a 2D second-order Taylor surface with residual below 1.5%. Two contributions to the principal-curvature splitting are decomposed. The oblique-projection component yields a free-space geometric astigmatic ratio determined solely by the arrival elevation, while ionospheric refraction acts as an anisotropic cylindrical lens, compressing the longitudinal curvature radius and expanding the transverse one. The phase mismatch between the simulated anisotropic wavefront and a single-radius isotropic reference grows quadratically with aperture extent. A closed-form expression for the aperture at which this mismatch reaches a half-power coherence-degradation threshold yields a scale of approximately 9 km for the present configuration. These results are presented as a configuration-specific case study whose analytic forms are intended to be portable.

Graphical abstract

Keywords

high-frequency(HF) skywave propagation / wavefront astigmatism / ionospheric ray tracing / deterministic phase-match aperture

1 Introduction

High-frequency skywave propagation utilizes the dynamic refractive properties of the Earth’s ionosphere to establish beyond-line-of-sight signal paths, serving as the critical physical carrier for transhorizon electronic reconnaissance, space-weather diagnostics, and passive emitter localization [1]. While single-site receiving architectures suffer from inherent resolution bounds and susceptibility to multi-path fading, deploying distributed passive ground arrays across baselines spanning several kilometers to tens of kilometers significantly enhances azimuthal resolution [2]. This distributed long-baseline paradigm capitalizes on the extensive spatial diversity of the incoming wavefield, fundamentally advancing the geolocation precision of distant uncooperative targets.

The lineage of far-field passive localization has traditionally relied on direction-of-arrival, time-difference-of-arrival, and frequency-difference-of-arrival estimation techniques [2]. These methods are anchored in the far-field plane-wave assumption, which posits a linear phase progression across the sensor array [3]. Its validity requires the source-to-array distance to exceed the Rayleigh distance 2D^2/\lambda [4]. For a typical HF skywave scenario (\lambda =30 m, distributed baseline D=10 km), the Rayleigh distance expands to approximately 6700 km, far exceeding the typical ionospheric slant range of several hundred kilometers. The distributed array therefore operates deeply within the near-field [5], where the plane-wave assumption and linear phase-steering techniques become inadequate.

To overcome the parameter coupling and data association ambiguities inherent in traditional two-step localization frameworks, direct position determination and passive synthetic aperture methodologies have seen aggressive development over the past two decades. Initiated by Weiss, direct position determination bypasses intermediate parameter extraction to construct a spatial cost function directly from baseband signals [6], later expanding to incorporate expectation-maximization algorithms for non-cooperative multi-source environments [7] and coprime array geometries for enhanced degrees of freedom [8]. In parallel, passive synthetic aperture techniques have demonstrated order-of-magnitude improvements in localization accuracy by leveraging relative platform motion [9], with recent advances effectively mitigating range-migration effects for linear frequency modulated radar sources [10]. Despite these robust algorithmic achievements, the ultimate geolocation ceiling of both direct position determination and passive synthetic aperture architectures remains strictly bottlenecked by the physical fidelity of their underlying wavefront phase models.

Providing high-fidelity phase histories for these distributed arrays requires advanced ionospheric ray tracing and propagation modeling. Geometric-optics ray tracing, originally formulated through the Haselgrove equations [11], has matured into comprehensive computational engines such as the IONORT software, which has recently been updated to support 3D global electron density grids [12,13]. Concurrently, to solve the specific boundary-value problem of emitter-to-receiver targeting, variational point-to-point alternatives have been developed [14,15]. Recently, generalized-force methods have advanced these variational approaches by computing the Hessian matrix of the phase-distance functional within a 3D ray-configuration space, enabling the precise mathematical classification of high and low rays via eigenvalue signatures [16,17]. It is crucial to distinguish, however, that these Hessian-based extensions operate exclusively within the configuration space for mode identification, rather than evaluating the 2D spatial phase surface at a physical receiver plane to extract principal wavefront curvature radii.

The literature addressing near-field arrays and high-frequency wavefront curvature exhibits a distinct gap between deterministic spatial modeling and statistical coherence analysis. Array signal processing conventionally assumes an isotropic quadratic phase structure governed by a single curvature radius [18,19], while the high-frequency propagation literature formalizes stochastic coherence lengths due to random refractive-index fluctuations [20,21]. However, classical optics demonstrates that oblique incidence upon a curved refractive boundary inevitably generates astigmatic aberrations, splitting the focal length into distinct sagittal and tangential components [22]. Neither the deterministic array models nor the stochastic scintillation theories treat the oblique skywave wavefront as a deterministic 2D anisotropic curvature-tensor problem. Specifically, the joint treatment of the geometric structural astigmatism induced by oblique projection and the anisotropic refractive lens modulation caused by along-path electron-density gradients remains unaddressed – this is the primary focus of the present study.

This paper formalizes the oblique high-frequency(HF) sky-wave wavefront as a deterministic 2D anisotropic curvature-tensor problem and quantifies its two underlying mechanisms on a representative single-hop oblique link. A virtual-source unfolded geometry is constructed from the IONORT-derived group path, and the receiver-plane phase surface is reconstructed from a dense ray bundle propagated through an IRI-2020 electron-density volume. Hessian-eigenvalue extraction applied separately to a free-space baseline and to the full ionosphere-loaded simulation isolates two distinct contributions: a closed-form geometric astigmatism arising purely from oblique projection, and an additional ionospheric lens modulation due to electron-density gradients. A closed-form Fresnel-integral coherence envelope is derived, yielding a deterministic phase-match aperture L_{\mathrm{c}\mathrm{o}\mathrm{h}} that bounds the spatial extent over which a single-radius wavefront model remains consistent with the simulated astigmatic wavefront.

2 System Geometry and Intrinsic Astigmatism

2.1 Oblique Skywave Link Geometry

A right-handed coordinate frame is established with its origin at the receiver (RX) reference point on the ground. The longitudinal axis \hat{\mathit{\xi }} is aligned with the great-circle bearing from RX toward the transmitter (TX); the transverse axis \hat{\mathit{\eta }} lies in the horizontal plane orthogonal to \hat{\mathit{\xi }}; and the vertical axis \hat{\mathit{\zeta }} is the local upward normal. TX is located at {39.9}^{\circ }N, {116.4}^{\circ }E (Beijing) and the RX reference point is placed at the landing centroid of the ray bundle, approximately {36.2}^{\circ }N, {115.5}^{\circ }E, corresponding to a great-circle separation D\approx 419 km (computed using standard spherical geometry).

For the wavefront-curvature analysis, the curved trajectory is represented by an unfolded straight-ray equivalent using a virtual source P whose position is fixed by the IONORT-derived group path P^{\prime } and the arrival elevation \alpha. Equivalently, P is constructed as the mirror image of the transmitter through the receiver tangent plane, displaced by the slant range R_0=D/\cos \alpha along the ray-arrival direction. This construction differs from the standard Breit–Tuve mid-point virtual height in that P is anchored to the receiver-side endpoint of the unfolded ray rather than to the transmitter-receiver mid-point, and serves purely as a geometric phase-source location. Denoting the one-way group path by P^{\prime }=c{t}_g, the corresponding virtual point P=({\xi }_P,0,{\zeta }_P) lies on the great-circle plane, where {\xi }_P is the along-track horizontal offset and {\zeta }_P is the vertical depth below the receiver tangent plane. For the reference link, {\zeta }_P=D\tan \alpha \approx 253 km. The slant range from P to the origin is

(1)R_0=({\xi }_P^2+{\zeta }_P^2{)}^{1/2}

(2)\sin \alpha ={\zeta }_P/R_0,\cos \alpha ={\xi }_P/R_0

The azimuth satisfies \theta =0 by the on-great-circle symmetry of the nominal path. The slant range R_0=D/\cos \alpha is the geometric straight-line distance from P to the receiver origin. It differs from the group path P^{\prime }=c{t}_g (the integrated optical path along the curved ray): ionospheric refraction (n_g>1) gives P^{\prime }>R_0, with the dispersive difference analysed in Section 3.3.

2.2 Isotropic-Spherical Reference Wavefront

Throughout this paper, the local horizontal coordinates are aligned with the longitudinal–transverse frame defined in Section 2.1: the x-axis coincides with \hat{\mathit{\xi }} (along the great-circle bearing toward the transmitter) and the y-axis coincides with \hat{\mathit{\eta }} (orthogonal to \hat{\mathit{\xi }} in the horizontal plane).

The monochromatic phase observed at an arbitrary horizontal offset \mathit{r}=(x,y,0{)}^{\mathsf{T}} relative to the origin, under the classical assumption of a locally spherical wavefront emanating from P, is

(3){\varphi }_{\mathrm{i}\mathrm{s}\mathrm{o}}(\mathit{r})=\frac{2\pi }{\lambda }(R(\mathit{r})-R_0)

where R(\mathit{r})=|\mathit{P}-\mathit{r}| is the exact slant distance. The canonical near-field expansion retains the second-order correction under a single scalar-radius (isotropic-spherical) assumption:

(4)R(\mathit{r})\approx R_0-\cos \alpha x+\frac{x^2{\sin }^2\alpha +y^2}{2R_0}

In (4), a single scalar radius R_0 governs the quadratic phase in both directions; the model admits no independent curvature parameter per axis. This single-parameter constraint is the structural limitation addressed in this paper.

2.3 Structural Astigmatism from Oblique Projection

The astigmatic structure of the geometrically propagated wavefront is revealed by evaluating the Hessian of R(\mathit{r}) at the origin. For \mathit{P}=({\xi }_P,0,{\zeta }_P), direct computation yields

(5){\frac{{\mathrm{\partial }}^{2}R}{\mathrm{\partial }{x}^2}|}_{\mathbf{0}}=\frac{{\sin }^2\alpha }{R_0},{\frac{{\mathrm{\partial }}^{2}R}{\mathrm{\partial }{y}^2}|}_{\mathbf{0}}=\frac{1}{R_0}

The effective principal curvature radii are therefore

(6)R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{g}\mathrm{e}\mathrm{o}}=\frac{R_0}{{\sin }^2\alpha },R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{g}\mathrm{e}\mathrm{o}}=R_0

With R_0=489 km and \alpha \approx {31.1}^{\circ } (Section 3.2), Eq. (6) give

(7)\begin{array}{rl}{R}_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{g}\mathrm{e}\mathrm{o}}& =1830\mathrm{k}\mathrm{m},R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{g}\mathrm{e}\mathrm{o}}=489\mathrm{k}\mathrm{m}\\ A{R}_{\mathrm{g}\mathrm{e}\mathrm{o}}& =\frac{R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{g}\mathrm{e}\mathrm{o}}}{R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{g}\mathrm{e}\mathrm{o}}}=3.74\end{array}

Two properties bear emphasis. First, R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{g}\mathrm{e}\mathrm{o}}\gg R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{g}\mathrm{e}\mathrm{o}}: the longitudinal curvature radius exceeds the transverse by a factor of 3.74, a direct consequence of oblique projection compressing the apparent longitudinal curvature () while leaving the transverse curvature unchanged. Second, {\mathrm{A}\mathrm{R}}_{\mathrm{g}\mathrm{e}\mathrm{o}}=1/{\sin }^2\alpha depends only on \alpha, which is fixed by the link geometry (ground range D and equivalent virtual height); it is independent of ionospheric state, frequency, or electron-density gradient. The closed-form 1/{\sin }^2\alpha prediction is validated against free-space Hessian fits in Section 3.3, and the dependence of the astigmatic ratio on arrival elevation — together with the ionosphere-loaded values — is presented for four benchmark links in Section 3.5.

3 Ionospheric Wavefront Modulation via Three-Dimensional Ray Tracing

3.1 Simulation Setup

The ionospheric electron density distribution is specified by the International Reference Ionosphere 2020 (IRI-2020) model [23] evaluated at epoch 2026-04-14T06:00 UTC, corresponding to daytime conditions over the link midpoint. A three-dimensional grid spanning {18}^{\circ }–{56}^{\circ }N latitude, {70}^{\circ }–{141}^{\circ }E longitude, and 1–600 km altitude is precomputed at 1^{\circ }\times 1^{\circ } horizontal and 2 km vertical resolution, and supplied as volumetric input to IONORT [12] (Fig. 1). Rays are traced at 10 MHz for the ordinary (O) polarization mode, with TX at {39.9}^{\circ }N, {116.4}^{\circ }E and RX near {36.2}^{\circ }N, {115.5}^{\circ }E, corresponding to a ground range D\approx 419 km and arrival elevation \alpha \approx {31.1}^{\circ }. The Hamiltonian ODE system is integrated with a relative error tolerance of 5\times {10}^{-4}, placing the typical F-layer step near 0.7 km and keeping the accumulated phase error below the level validated in Section 3.2. All rays flagged as success additionally satisfy a ground-range criterion D_{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}>100 km, which eliminates transmitter-collocated truncation pseudo-rays arising near the MUF watershed; across the four reference links the minimum landing range is 360 km, comfortably above this threshold. IONORT ray tracing indicates that the rays turn around at approximately 122 km altitude. As the wavefront analysis of Sections 2–4 focuses on the integrated 2D phase structure at the receiving plane, it remains independent of which specific ionospheric layer produces the reflection; the curvature extraction procedure applies to the actual simulated propagation regardless of turning height.

3.2 2D Wavefront Surface Reconstruction and Hessian Matrix Extraction

Wavefront reconstruction proceeds by launching a dense ray bundle over a 51\times 61 elevation–azimuth grid (elevation {26.88}^{\circ }–{36.88}^{\circ } at {0.2}^{\circ } step; azimuth {187.43}^{\circ }–{193.43}^{\circ } at {0.1}^{\circ } step). Of 3111 launched rays, 2748 satisfy the success criterion and terminate on the horizontal sampling plane within an irregular footprint of approximately 65\times 130 km centred on the L1 reference point. The accumulated propagation phase \mathrm{\Phi }(x_k,y_k) is recorded at each landing point and fitted to a 2D second-order Taylor expansion:

(8)\mathrm{\Phi }(x,y)={\mathrm{\Phi }}_0+g_{\xi }x+g_{\eta }y+{\displaystyle \frac{1}{2}}{H}_{xx}{x}^2+H_{xy}xy+{\displaystyle \frac{1}{2}}{H}_{yy}{y}^2

with six coefficients determined by least-squares regression over the 2748 scattered landing positions. The symmetric phase-Hessian

(9)\mathbf{H}=\left(\begin{array}{cc}{H}_{xx}& H_{xy}\\ H_{xy}& H_{yy}\end{array}\right)

is diagonalized to align the principal axes with \hat{\mathit{\xi }} and \hat{\mathit{\eta }}. The two principal curvature radii are

(10)R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}}=\frac{2\pi /\lambda }{H_{xx}},R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}=\frac{2\pi /\lambda }{H_{yy}}.

The maximum pointwise residual across all 2748 landing points is less than 1.5%, confirming that the second-order Taylor model captures the dominant curvature structure with negligible truncation error.

The 2748 landing points form a non-uniform spatial density (denser near the centroid). A bootstrap test (1000 draws, 80% subsampling) gives 95% confidence interval half-widths of 0.25% and 0.66% for the longitudinal and transverse curvature radii, with {\mathrm{A}\mathrm{R}}_{\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}} varying under 1%. The fit is therefore robust to the non-uniform sampling, owing to the smoothness of the quadratic phase surface and the heavily overdetermined regression (six parameters on 2748 points).

The Hessian fit is performed only on the 2748successfully landed rays. The reconstructed phase surface and fit overlay are shown in Fig. 2.

3.3 Anisotropic Ionospheric Lens Effect

The Hessian extraction is applied separately to a geometry-only baseline (ionosphere disabled) and to the full IONORT simulation at 10 MHz, O mode. For the geometry-only baseline, the same Taylor-fit pipeline recovers R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{g}\mathrm{e}\mathrm{o}}^{\mathrm{f}\mathrm{i}\mathrm{t}}=1750 km and R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{g}\mathrm{e}\mathrm{o}}^{\mathrm{f}\mathrm{i}\mathrm{t}}=475 km, deviating by 4.4\mathrm{\%} and 2.9\mathrm{\%} respectively from the analytic predictions of (6) (1830 km and 489 km). The fitted astigmatic ratio {\mathrm{A}\mathrm{R}}_{\mathrm{g}\mathrm{e}\mathrm{o}}^{\mathrm{f}\mathrm{i}\mathrm{t}}=3.69 differs from the theoretical 3.74 by 1.5\mathrm{\%}. These small discrepancies reflect the finite, irregular sampling of the landing footprint and the tangent-plane coordinate projection used identically for the ionosphere-loaded case, confirming that the pipeline recovers the closed-form geometry to within a few percent.

Ionospheric refraction modifies both principal radii in opposite senses. The longitudinal radius is compressed to R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}=1388 km (a 20.7% reduction), reflecting convergence driven by the along-path electron-density gradient through the reflection region. The transverse radius is expanded to R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}=537 km (a 13.0% increase), reflecting weak transverse divergence from the cross-path gradient. The ionospheric astigmatic ratio is {\mathrm{A}\mathrm{R}}_{\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}=2.59.

Along the longitudinal (great-circle) direction, the along-path electron-density gradient produces gradient-index refraction that bends rays convergently toward the ground plane, equivalent to a converging cylindrical lens whose focal power is set by the along-path N_e gradient. Along the transverse direction, the cross-path gradient is both weaker and directed outward, producing mild divergence equivalent to a weakly diverging lens. The combined action reduces AR from 3.74 to 2.59: the ionosphere partially mitigates the intrinsic geometric astigmatism, but the structural factor 1/{\sin }^2\alpha persists regardless of ionospheric state.

Dispersive medium response and positioning implications. The IRI-2020+IONORT pipeline inherently models the ionosphere as a dispersive magnetised plasma. In the HF band, the Appleton–Hartree relation [24] gives and n_g=\mathrm{d}(\omega n_p)/\mathrm{d}\omega >1. Integrated along the ray trajectory \mathrm{\Gamma },

(11)P^{\prime }={\int }_{\mathrm{\Gamma }}{n}_g\mathrm{d}s>{\int }_{\mathrm{\Gamma }}\mathrm{d}s>{\int }_{\mathrm{\Gamma }}{n}_p\mathrm{d}s={\mathrm{\Phi }}_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{h}}

with strict inequalities holding if and only if \mathrm{\partial }{n}_p/\mathrm{\partial }\omega \ne 0. Across the 2748 successfully landed rays, ⟨P^{\prime }-{\mathrm{\Phi }}_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{h}}⟩\approx 24.33 km (peak 31.13 km). Analysis of the most central ray (Section 3.2, selection criterion: smallest centroid distance) shows the cumulative dispersive difference {\int }_{}^{}(n_g-n_p)\mathrm{d}s=23.9 km, corresponding to a group delay of 0.080 ms. This localization has practical implications for HF positioning: the group path exceeding the phase path by tens of kilometers implies that time-of-arrival and phase-based ranging yield systematically different apparent distances, and the concentration of dispersive delay in a \sim 15 km turning-layer slab suggests that correction models need only capture the electron-density structure of this compact altitude range rather than the full profile. The anisotropic curvature modulation extracted above therefore rests upon an explicitly dispersive refractive-index field.

3.4 Cross-Link Generalisation

The full pipeline is re-applied to three additional oblique links spanning \alpha \in [{6.7}^{\circ },{31.1}^{\circ }] and f\in [10,18] MHz, with ground distances D from 419 km to 1414 km, sharing the common epoch and O-mode single-hop regime. The L3 carrier f=18 MHz lies well above f_{\mathrm{o}}F2\approx 9 MHz but safely within the oblique MUF: at \alpha ={6.74}^{\circ } the secant-law obliquity factor \sec {\varphi }_i\approx 3.4 yields a single-hop MUF on the order of 30 MHz, consistent with the 1281 successful rays returned by IONORT. At ground ranges exceeding approximately 1000 km (L2, L3), the tangent-plane Hessian includes a parasitic quadratic contribution from Earth curvature, introducing a systematic bias in the extracted R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}} not de-embedded here. This bias arises from the departure of the local tangent plane from the spherical Earth: for R_E\approx 6371 km, the tangent-plane approximation introduces a quadratic term \sim D^2/(2R_E) in the apparent longitudinal curvature, negligible at L1 () but reaching several percent at L3. De-embedding this contribution requires a spherical-coordinate Hessian formulation, left for future work. The benchmark is summarised in Tab. 1 and visualised in Fig. 3 and Fig. 4.

Three observations emerge. First, the geometric astigmatic ratio follows 1/{\sin }^2\alpha across the full elevation sweep, with Taylor-fit values within 14% of the closed-form limit even at the most grazing link (\alpha ={6.74}^{\circ }, {\mathrm{A}\mathrm{R}}_{\mathrm{g}\mathrm{e}\mathrm{o}}^{\mathrm{f}\mathrm{i}\mathrm{t}}=59.9 versus 1/{\sin }^2\alpha =72.6). Second, ionospheric refraction consistently reduces AR, with the reduction growing monotonically at lower elevation: from 30% at \alpha \approx {31}^{\circ } (3.74 \to 2.59) to 48% at \alpha \approx {6.7}^{\circ } (59.9 \to 31.3), reflecting the longer along-path residence within the stratified electron-density gradient at grazing incidence. Third, even at the lowest elevation {\mathrm{A}\mathrm{R}}_{\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}} remains bounded away from unity by more than an order of magnitude, confirming that 2D astigmatism rather than isotropic sphericity is the dominant wavefront geometry.

4 Phase Mismatch Budget and Deterministic Phase-Match Aperture

4.1 Residual Phase Mismatch Against an Isotropic-Spherical Reference

The true wavefront phase is encoded by \mathit{H}, whose principal curvatures are 1/R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}} and 1/R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}. A radio-measurement model adopting the isotropic-spherical reference (4) with scalar radius R_0 incurs a residual phase mismatch. Defining the degenerate reference Hessian {\mathit{H}}_{\mathrm{i}\mathrm{s}\mathrm{o}}=(2\pi /\lambda R_0)\mathit{I}, the mismatch is

(12)\begin{array}{rl}\mathrm{\Delta }\psi (\mathit{r})=& \frac{1}{2}{\mathit{r}}^{\mathsf{T}}({\mathit{H}}_{\mathrm{i}\mathrm{s}\mathrm{o}}-\mathit{H})\mathit{r}=\\ & \frac{\pi }{\lambda }[x^2\left(\frac{1}{R_0}-\frac{1}{R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}}\right)+\\ & y^2\left(\frac{1}{R_0}-\frac{1}{R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}}\right)]\end{array}

where R_0=489 km is adopted as the reference radius. Along the longitudinal cut (y=0),

(13)\mathrm{\Delta }{\psi }_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}}(x)=\frac{\pi }{\lambda }{x}^2\left(\frac{1}{R_0}-\frac{1}{R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}}\right)

and along the transverse cut (x=0),

(14)\mathrm{\Delta }{\psi }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}(y)=\frac{\pi }{\lambda }{y}^2\left(\frac{1}{R_0}-\frac{1}{R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}}\right)

At equal offset,

(15){\frac{\mathrm{\Delta }{\psi }_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}}}{\mathrm{\Delta }{\psi }_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}}|}_{|x|=|y|}=\frac{1/489-1/1388}{1/489-1/537}=7.3

reflecting an 899 km longitudinal underestimate of the curvature radius by the isotropic reference versus only a 48 km transverse underestimate. The mismatch profiles are shown in Fig. 5.

4.2 Derivation of the Deterministic Phase-Match Aperture L_{\mathrm{c}\mathrm{o}\mathrm{h}}

The spatial scale over which the isotropic reference remains consistent with the true astigmatic wavefront is bounded by accumulation of the quadratic mismatch. Consider the longitudinal direction (the dominant mismatch); the transverse case follows by substituting R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}.

Define the longitudinal curvature mismatch \mathrm{\Delta }\triangleq 1/R_0-1/R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}. Operating in a phase-dominant regime – the spatial amplitude variations of the ray-ensemble wavefront across the local sampling window are small compared with the quadratic phase mismatch at the aperture scales considered here – the normalised aperture-integrated coherence amplitude across [-L/2,L/2] is

(16)G(L)=\frac{1}{L}\left|{\int }_{-L/2}^{L/2}{\mathrm{e}}^{\mathrm{j}(\pi /\lambda )\mathrm{\Delta }{x}^2}\mathrm{d}x\right|

Introducing the dimensionless half-aperture t(L)=(L/2)\sqrt{2|\mathrm{\Delta }|/\lambda }, (16) evaluates to a normalized Fresnel integral. The deterministic phase-match aperture L_{\mathrm{c}\mathrm{o}\mathrm{h}} is defined as the smallest sampling extent at which the isotropic reference departs from the true wavefront by half-power (\mathcal{L}=3 dB, G=1/\sqrt{2}). Evaluating the Fresnel integral yields the universal dimensionless threshold t_{3\mathrm{d}\mathrm{B}}\approx 1.32, giving

(17)L_{\mathrm{c}\mathrm{o}\mathrm{h}}=2t_{3\mathrm{d}\mathrm{B}}\sqrt{\frac{\lambda }{2|\mathrm{\Delta }|}}=1.864\sqrt{\frac{\lambda }{|\mathrm{\Delta }|}}

Substituting \lambda =30 m and |\mathrm{\Delta }|=|1/489-1/1388|{\mathrm{k}\mathrm{m}}^{-1}=1.324\times {10}^{-6}{\mathrm{m}}^{-1} yields

(18)L_{\mathrm{c}\mathrm{o}\mathrm{h}}=8.87\mathrm{k}\mathrm{m}

The transverse coherence length, obtained by substituting R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}=537 km, is L_{\mathrm{c}\mathrm{o}\mathrm{h},\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}=24.0 km. The factor-of-2.7 ratio between the two directional coherence lengths mirrors the factor-of-7.3 disparity in quadratic phase accumulation (15), inverted by the square-root in (17). These two scales define the anisotropic validity envelope of a single-radius wavefront description for the analysed link.

4.3 Cross-Link Consistency of the Longitudinal Aperture Bound

Re-applying (17) to the four-link benchmark yields the values in Tab. 1. Across the four links, L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\parallel } remains within 8.87–10.50 km (mean 9.40 km, spread \pm 9\mathrm{\%}) despite {\mathrm{A}\mathrm{R}}_{\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}} spanning a 12\times dynamic range. Inspection of (17) explains this: as elevation decreases, R_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}} grows super-linearly, so |\mathrm{\Delta }| asymptotically approaches 1/R_0 alone, and the \sqrt{\lambda /|\mathrm{\Delta }|} scaling absorbs the residual variation. The transverse bound L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\perp } exhibits a wider spread (23.1–30.5 km), reflecting weaker curvature modulation. Within the limits of this study, the four-link consistency suggests L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\parallel } may serve as an order-of-magnitude reference for sub-aperture sizing.

5 Implications for Distributed Long-Baseline HF Passive Localization

The deterministic phase-match aperture L_{\mathrm{c}\mathrm{o}\mathrm{h}} characterises the spatial scale over which a single-radius isotropic-spherical steering model remains internally consistent with the simulated 2D astigmatic wavefront. This section translates that scale into theoretical bounds relevant to distributed long-baseline HF receiver arrays. The estimates below are order-of-magnitude consequences of the present wavefront model, not absolute predictions of fielded system accuracy.

5.1 Local Direction-of-Arrival Disagreement

A distributed coherent aperture composed of sub-aperture nodes samples the arriving wavefront at different offsets. A sub-aperture at longitudinal offset L/2 from the origin samples a local phase gradient that, under the isotropic model, yields a local wavefront normal differing from the true astigmatic one. Differentiating (13) gives the local gradient, which maps to an elevation-DOA disagreement

(19)\delta \alpha (L)\approx \frac{|\mathrm{\Delta }|L}{2}

Evaluating at L=L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\parallel }=8.87 km gives \delta \alpha \approx {0.34}^{\circ }. The transverse counterpart, with R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s},\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}} and L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\perp }=24.0 km, gives \delta \beta \approx {0.12}^{\circ }. Below these bounds the disagreement decreases linearly with L; above them the quadratic phase residual exceeds the \pi /2 threshold and the linear-tilt approximation ceases to hold. Within a single sub-aperture, the same quadratic mismatch produces wavefront defocus (main-lobe broadening) rather than a DOA bias. Fig. 6 plots \delta \alpha (L) and \delta \beta (L) with the aperture bounds marked.

5.2 Ground-Localization Bias

For a single-station skywave geometry, the ground range D\approx {\zeta }_P/\tan \alpha and its sensitivity to elevation-DOA error is |\mathrm{\partial }D/\mathrm{\partial }\alpha |={\zeta }_P/{\sin }^2\alpha. Substituting {\zeta }_P\approx 253 km and \alpha \approx {31.1}^{\circ } yields |\mathrm{\partial }D/\mathrm{\partial }\alpha |\approx 947 km/rad. The resulting range-direction displacement at the longitudinal aperture bound is

(20)\delta D(L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\parallel })\sim 5.5\mathrm{k}\mathrm{m}

The transverse DOA disagreement maps to a cross-range displacement through the great-circle distance: \delta D_{\perp }\approx D\delta \beta, giving \delta D_{\perp }(L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\perp })\sim 1 km. These estimates are order-of-magnitude theoretical extrapolations consistent with the present wavefront model.

5.3 Indications for Distributed Array Design

The longitudinal aperture bound L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\parallel }\approx 8.87 km marks the scale beyond which an isotropic-spherical steering vector incurs more than -3 dB coherent-integration penalty against the true 2D astigmatic wavefront. Within a single coherent sub-aperture, longitudinal extents on the order of L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\parallel } (kilometers) may serve as a useful sizing reference; the transverse extent is bounded by L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\perp }\approx 24.0 km. Coherent processing across sub-apertures separated by baselines exceeding L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\parallel } would require steering vectors that explicitly encode the 2D anisotropic near-field phase characterised in Section 3 – parameterised jointly by longitudinal and transverse principal-curvature radii – together with accommodation of the geometric Jacobian |\mathrm{\partial }D/\mathrm{\partial }\alpha |. It is important to note that L_{\mathrm{c}\mathrm{o}\mathrm{h}} derived here is a deterministic bound based on the smooth IRI-2020 mean ionosphere; in practice, traveling ionospheric disturbances (TIDs) and small-scale scintillation will further degrade spatial coherence at scales smaller than L_{\mathrm{c}\mathrm{o}\mathrm{h}}, so the present bound should be interpreted as the model's upper limit of systematic wavefront error rather than the coherence floor of the physical channel. It should also be noted that the wavefront analysis developed here characterises the arriving field with respect to the virtual source P; relating this virtual-source location to the physical transmitter coordinates involves the inverse ionospheric ray-tracing problem and is deferred to subsequent work. A complete treatment of such anisotropic steering vectors and associated algorithms lies outside the descriptive scope of this paper and is left as future work.

6 Conclusion

This paper reported a simulation study of the two-dimensional astigmatic phase structure of HF skywave wavefronts on a short single-hop oblique link (transmitter at Beijing, ground range D\approx 419 km, 10 MHz, O mode, IRI-2020 epoch 2026-04-14T06:00 UTC), separating geometric and ionospheric contributions through IONORT 3D ray tracing and 2D Hessian extraction with fit residual below 1.5%.

Within the scope of this study, three findings are reported. First, pure oblique geometry produces an astigmatic ratio {\mathrm{A}\mathrm{R}}_{\mathrm{g}\mathrm{e}\mathrm{o}}=3.74, governed by the arrival elevation through {\mathrm{A}\mathrm{R}}_{\mathrm{g}\mathrm{e}\mathrm{o}}=1/{\sin }^2\alpha, independent of ionospheric state. Second, ionospheric refraction is consistent with an anisotropic cylindrical lens: longitudinal convergence of 20.7% and transverse divergence of 13.0% reshape the ratio to {\mathrm{A}\mathrm{R}}_{\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}}=2.59. The associated phase mismatch against an isotropic-spherical reference exceeds the half-power threshold at a deterministic phase-match aperture L_{\mathrm{c}\mathrm{o}\mathrm{h}}=8.87 km. Extending the analysis to four oblique links (\alpha \in [{6.7}^{\circ },{31.1}^{\circ }], f\in [10,18] MHz) shows {\mathrm{A}\mathrm{R}}_{\mathrm{g}\mathrm{e}\mathrm{o}} tracking the 1/{\sin }^2\alpha prediction within 14% and L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\parallel } remaining within 9.4\pm 0.9 km. Third, the present simulation suggests that this \sim 10 km longitudinal aperture scale may serve as a useful sub-aperture sizing reference for distributed long-baseline HF receiver arrays: coherent integration within a single sub-aperture of longitudinal extent \lesssim L_{\mathrm{c}\mathrm{o}\mathrm{h}}^{\parallel } can proceed under an isotropic-spherical steering model without exceeding the -3 dB coherence-degradation threshold.

Extension to other diurnal/seasonal conditions, geographic regions, geomagnetic disturbance levels, and propagation modes is left as future work.

References

[1] S. J. Anderson . Societal applications of hf skywave radar. Remote Sensing, 2022, 14(24): 6287

[2] G. Wu , F. Guo , M. Zhang . Direct position determination: An overview. Journal of Radars, 2020, 9(6): 998–1013

[3] H. Krim , M. Viberg . Two decades of array signal processing research: The parametric approach. IEEE Signal Processing Magazine, 1996, 13(4): 67–94

[4] A. M. Elbir , T. E. Tuncer . Far-field doa estimation and near-field localization for multipath signals. Radio Science, 2014, 49(9): 765–776

[5] V. E. Gherm , N. N. Zernov . Fresnel filtering in hf ionospheric reflection channel. Radio Science, 1995, 30(1): 127–134

[6] A. J. Weiss . Direct position determination of narrowband radio frequency transmitters. IEEE Signal Processing Letters, 2004, 11(5): 513–516

[7] E. Tzoreff , A. J. Weiss . Expectation-maximization algorithm for direct position determination. Signal Processing, 2017, 133: 32–39

[8] G. Wu , Y. Zhang , W. Zhang , Y. Pei , M. Zhang , F. Guo . Coprime array based direct position determination of signals with single moving observation. Journal of Radars, 2022, 11(4): 692–704

[9] A. Li , H. Huan , R. Tao , L. Zhang . Passive synthetic aperture high-precision radiation source location by single satellite. IEEE Geoscience and Remote Sensing Letters, 2022, 19: 1–5

[10] D. Li , D. Wang , H. Huan . LFM radar source passive localization algorithm based on range migration. Journal of Beijing Institute of Technology, 2024, 33(2): 130–140

[11] J. Haselgrove, “Ray theory and a new method for ray tracing,” in Proceedings of the Cambridge Conference on the Physics of the Ionosphere, pp. 355–364, 1955.

[12] A. Azzarone , C. Bianchi , M. Pezzopane , M. Pietrella , C. Scotto , A. Settimil . Ionort: A windows software tool to calculate the hf ray tracing in the ionosphere. Computers & Geosciences, 2012, 42: 57–63

[13] M. Pietrella , M. Pezzopane , A. Pignatelli , A. Pignalberi , A. Settimi . An updating of the ionort tool to perform a high-frequency ionospheric ray tracing. Remote Sensing, 2023, 15: 21

[14] C. J. Coleman . A ray tracing formulation and its application to some problems in over-the-horizon radar. Radio Science, 1998, 33(4): 1187–1197

[15] C. J. Coleman . Point-to-point ionospheric ray tracing by a direct variational method. Radio Science, 2011, 46(5): RS5016

[16] I. A. Nosikov , M. V. Klimenko , G. A. Zhbankov , A. V. Podlesnyi , V. A. Ivanova , P. F. Bessarab . Generalized force approach to point-to-point ionospheric ray tracing and systematic identification of high and low rays. IEEE Transactions on Antennas and Propagation, 2020, 68(1): 455–467

[17] I. A. Nosikov , M. V. Klimenko , A. M. Padokhin , V. V. Nosikova , P. F. Bessarab . Generalized force method for point-to-point ray tracing in anisotropic ionosphere: Implementation and applications to Nequick2 and IGRF13 models. Radio Science, 2025, 60(3): e2024RS008092

[18] B. Friedlander . Localization of signals in the near-field of an antenna array. IEEE Transactions on Signal Processing, 2019, 67(15): 3885–3893

[19] F. Lu , H. Zhao , X. Zhao , X. Wang , A. Saleem , G. Zheng . Investigation of near-field source localization using uniform rectangular array. Electronics, 2022, 11(12): 1916

[20] K. C. Yeh , C. H. Liu , M. Y. Youakim . A theoretical study of the ionospheric scintillation behavior caused by multiple scattering. Radio Science, 1975, 10(1): 97–106

[21] H. G. Booker . The role of acoustic gravity waves in the generation of spread-f and ionospheric scintillation. Journal of Atmospheric and Terrestrial Physics, 1979, 41(5): 501–515

[22] M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, 7 th ed. Cambridge: Cambridge University Press, 1999.

[23] D. Bilitza , M. Pezzopane , V. Truhlik , D. Altadill , B. W. Reinisch , A. Pignalberi . The International Reference Ionosphere model: A review and description of an improved IRI-2020 version. Journal of Space Weather and Space Climate, 2022, 12: 25

[24] K. G. Budden, The Propagation of Radio Waves: The Theory of Radio Waves of Low Power in the Ionosphere and Magnetosphere. Cambridge: Cambridge University Press, 1985.

Rights & permissions

The associate editor coordinating the review of this manuscript was Dr. Junhua Yang