Introduction
Radio-frequency neutrino detection in cold polar ice
Over the last three decades, aided by significant advances in digital technology, detection of coherent radio-frequency (RF) Cherenkov-like radiation (Askar’yan, Reference Askar’yan1961, Reference Askar’yan1965), or RF echoes from neutrino-induced cascades (Prohira and others, Reference Prohira2021) in dense media has emerged as the most cost-effective approach to measure Ultra-High Energy Neutrinos (UHEN). Such UHEN detectors need to monitor 
$\mathcal{O}({10}\,\mathrm{km}^{3})$ volumes to be sensitive to the extremely low UHEN flux at neutrino energies 
$\sim 100\,\mathrm{PeV}$ and above. Abundant, high-purity polar ice provides a nearly ideal dense target medium, featuring 
$\mathcal{O}(1\,\mathrm{km})$ RF attenuation lengths (Robin and others, Reference Robin, Evans and Bailey1969; Barwick and others, Reference Barwick, Besson, Gorham and Saltzberg2005; Aguilar and others, Reference Aguilar2022), and typically moderate anthropogenic noise environments. Figure 1 illustrates the neutrino detection scheme; further reviews of UHEN detection using RF antennas can be found in the literature (Huege and Besson, Reference Huege and Besson2017; Schröder, Reference Schröder2017).
Schematic illustrating UHEN signal detection. An incident neutrino collides with an ice molecule at a ‘vertex’, producing Askaryan radiation (consisting of radially outward polarized propagating electric fields) geometrically confined to the surface of an expanding Cherenkov cone. As a result of the variable RI in the firn (deep purple), two rays (direct [‘D’; solid lines] and refracted/reflected [‘R’; dashed lines]) reach each of the Surface and also Deep receiver antennas. Purple lines show the highest-amplitude, ‘on-Cone’ signal and red lines show weaker, ‘off-Cone’ signal.
RF ice permittivity and the impact on UHEN experiments
Accurate and precise characterization of the RF response of glacial firn is essential for UHEN experiments seeking observation of neutrinos using radio techniques, as the dielectric permittivity in the hundreds of MHz frequency regime directly impacts the active volume over which neutrino interactions can be measured (‘effective volume’) and also the ability to point an interacting neutrino back into the sky (‘neutrino astronomy’). One parameter of particular interest is the depth-dependent refractive index profile (‘RIP’, or 
$n(z)$), which determines the speed of light in different media and therefore dictates signal trajectories. For Radio Neutrino Observatory-Greenland (RNO-G), with antennas deployed 1–20 m apart vertically, and 1.25 km laterally between receiver stations, 1% precision in local RIP measurements is required in order that 
$n(z)$ uncertainties not dominate neutrino-finding and reconstruction errors.
Refractive index (RI) modeling
The dielectric response of a given medium to an electromagnetic excitation of angular frequency 
$\omega$ is often formulated in the Drude–Lorentz atomic-resonator model (Drude, Reference Drude1900; Lorentz, Reference Lorentz1916), for which the equation of motion is calculated assuming an incident vector Electric field 
${\vec E}$ driving an atomic electron, modulo damping (via interatomic effects and giving rise to extinction) and resonance effects. Ice lattice interactions produce an absorptive resonance pole in the complex pure ice permittivity at far-infrared frequencies (4.8 and 14.48 THz). At RFs and below, and extending into the hundreds of kHz regime, ice response is characterized by temperature-dependent Debye relaxation, corresponding to small absorptive losses and therefore long attenuation lengths.
The Fujita et al. compilation of extant complex dielectric data (Fujita and others, Reference Fujita, Matsuoka, Ishida, Matsuoka and Mae2000) extensively summarizes both lab and in situ measurements of the complex ice permittivity, indicating a 
$\sim$0.5% systematic uncertainty in the RI, with a 0.2% variation over the temperature range relevant for UHEN experiments. Consistent with expectations from the atomic resonator model, the dependence of the complex permittivity on wavelength over the RF band is mild. In their review, the authors specifically consider the physical attributes of natural ice (density, air bubble infusion, impurity [particularly conducting acids and free 
$H^+$ ions] and temperature) affecting the refractive indices, finding density to be the dominant factor and supporting an assumed linear dependence of RI on density. Their review of the literature indicates a secondary role played by temperature; translating their compiled data into RI temperature dependence, their data imply 
$n(T)$=
$n(T_0$)+
$\alpha(T-T_0)$, with 
$n(T_0=190$K)=1.77 and 
$\alpha=8.3\times 10^{-5}$ over our relevant temperature range. There was no conclusive frequency dependence (dispersion) indicated by current experimental results; in fact, strong limits have been placed at both South Pole and also Summit Station, ruling out dispersion as an important contributor to UHEN experimental systematics (Kravchenko and others, Reference Kravchenko, Besson, Ramos and Remmers2011; Aguilar and others, Reference Aguilar2023).
RI measurements to date
The real part of the permittivity (
$\varepsilon_r$) can be directly quantified in the lab via dielectric profiling of ice cores, the small lateral size of the cores being profiled, relative to radio wavelengths, notwithstanding.
Gravimetric density measurements are also straightforward, and can be translated to permittivity assuming a one-parameter dependence (Schytt, Reference Schytt1958; Kovacs and others, Reference Kovacs, Gow and Morey1995):
\begin{align}
\varepsilon_r = \Bigg[1+0.845\rho\Big(\frac{\mathrm{g}}{\mathrm{cm}^{3}}\Big)\Bigg]^2;~~ n=\sqrt{\varepsilon_r}.
\end{align}Neutron probe monitor (NPM) measurements (Morris and Cooper, Reference Morris and Cooper2003; Hawley and Morris, Reference Hawley and Morris2006), made in an ice borehole for which the hole diameter has already been measured to a precision of 400 microns can, in principle, measure 2% density (statistical) variations on a vertical length scale of 1 cm if neutron back-scatter can be measured at each sampling point for 64 s (the hole diameter must be known in order to precisely calibrate out residual back-scatter unassociated with the ice volume itself).
The bulk index-of-refraction can also be derived in situ by the UHEN experiments using transit time (TT) measurements of polarized impulsive (ns-duration) signals propagating over a macroscopic distance d, for which TT=d
$n/c_0$. This technique was recently applied to surface-based radar-echo sounding (RES) data taken at Summit Station (Welling and others, Reference Welling2024), for which the time difference between successive internal layer echoes could be determined to nanosecond-scale accuracy. Correlating the echo times with decimeter-scale tabulations of internal layer depths from archived ice cores in the vicinity of Summit Station provided the most precise determination, to date, of the bulk polar ice RI (
$n$=1.778
$\pm$0.006) at depths below several hundred meters (Welling and others, Reference Welling2024).
RNO-G RIP parameterization and current limitations
Translations of density into 
$n(z)$ are one important component of the RIP used by the RNO-G experiment (Windischhofer, Reference Windischhofer2024). Performing a fit over all the available techniques and data yields the RIP favored by the RNO-G experiment (Windischhofer, Reference Windischhofer2024; Agarwal and others, Reference Agarwal2025). One major drawback of these fits, however, is that they yield a single (smooth) 
$n(z)$ parameterization for the entire 35 
${\mathrm{km}}^2$ RNO-G footprint. However, significant (5–10%) local variations in 
$n(z)$ have already been noted at sites less than 1 km away from each other (Morris and Cooper, Reference Morris and Cooper2003; Hawley and Morris, Reference Hawley and Morris2006) near Summit.
One-parameter 
$n(z)$ formulations are also insensitive to the (known) tensorial characteristics of the dielectric permittivity. Several experiments have demonstrated polarization-dependent anisotropies in RF propagation depending on the orientation of the electric field vector relative to the ice crystal orientation fabric (COF). As is well-known at optical frequencies, and also now well-measured at RFs (Kravchenko and others, Reference Kravchenko, Besson, Ramos and Remmers2011; Allison and others, Reference Allison2019; Besson and others, Reference Besson, Kravchenko and Nivedita2023), birefringence results in a wavespeed asymmetry for polarizations parallel vs perpendicular to the COF. Fujita et al. predict a 0.5% asymmetry, although direct in situ experiments have only measured (at most) half that value, perhaps as a result of incomplete alignment of the COF.
Additionally, (known) density fluctuations can disrupt the coherence between the frequency components that comprise the sought-after otherwise-coherent Cherenkov radiation signal electric field vector, for either neutrinos observed via Refracted signals with trajectories sampling the near-surface ice, or down-coming cosmic-ray-induced showers with steep inclination angles. Recently, the Askaryan Radio Array claimed the first-ever observation of such down-coming cosmic ray signals (Alden and others, Reference Alden2025). Numerically, we can estimate the magnitude of decoherence from density fluctuations, knowing that the number of particles contributing to the coherent Cherenkov signal is approximately the shower energy, in GeV, divided by 4. A shower having energy 
$10^{17}$ eV, or 
$10^8$ GeV, would therefore be reduced in amplitude by a factor of 10,000 in the extreme case of full coherence loss, underscoring the importance of mapping any density fluctuations on a length scale comparable to the smallest measurable wavelengths (20 cm) for the UHEN experiments.
Overdense melt layers within the firn, for example, can vary in thickness from centimeter to meter scale, with local density variations that can approach that of ice, and therefore produce significant coherence loss (Nghiem and others, Reference Nghiem2012). Quantifying the magnitude of density fluctuations (which are empirically observed to decrease in magnitude with increasing depth) can (hopefully) allow us to determine the impact on UHEN and UHECR astronomy (Determining the numerical implications of density fluctuations on measurements of Askaryan radiation from impacting air shower cores is a primary goal of an upcoming Summit Station field season).
Glacial firn ice properties
The antennas of the RNO-G experiment are deployed at maximum depths of 100 m, within the firn layer of the ice sheet. Firn properties change with depth, as the ice gradually compacts with increasing overburden. Unlike the Earth’s atmosphere, which readily lends itself to a simple exponential density profile for a fluid compressing under its own weight, compaction of accumulated snow proceeds through multiple density transitions, often demarcated by two separate inflection points (Herron and Langway, Reference Herron and Langway1980; Salamatin and others, Reference Salamatin, Lipenkov and Duval1997; Stevens and others, Reference Stevens2020). At a density of 550 
${\mathrm{km}\ \mathrm{m}^{-3}}$, corresponding to a depth of 12 m at Summit Station, snow (depending on moisture content) crystals transform to ice characterized by an abundance of grains and a density intermediate between snow and glacial ice. At a density of 830 
${\mathrm{km}\ \mathrm{m}^{-3}}$, ‘close-off’ occurs and bubbles can no longer mix with surface air. At this depth (
$\sim$75–80 m at Summit), firn ice transforms to ‘bubbly ice’, beyond which the ice structure is relatively constant, modulo a diminishing fraction of air inclusions. Over the upper 100 m, the density profile thus evolves from 
$\sim$45% of the asymptotic deep ice density (917 kg 
${\mathrm{m}^{-3}}$) to nearly 90–95% of the asymptotic value. This density change corresponds to an RI range of 
$\sim$1.3–1.75 at Summit; averaged over the entire Greenland Ice Sheet, the average surface (upper 10 cm) ice density has been measured to be 315
$\pm$44 kg 
${\mathrm{m}^{-3}}$ (Fausto and others, Reference Fausto2018), corresponding to a surface RI value 
$\textit{n}_{\text{surface}}(z)=1.27$.
Thus far, numerous models have been advanced to parameterize firn densification relative to the surface; each model should ideally include a quantitative description of multiple factors contributing to the observed 
$\rho(z)$ density profile, including local moisture content, local snowfall rate, temperature profile and heat diffusion through the ice, diffusion of both water and above-surface air, grain growth and evolution of the local ice fabric, etc. The large number of input priors underscores the importance of obtaining in situ density data (Vandecrux and others, Reference Vandecrux2020). More extensive information regarding ice, firn and scattering and densification processes can be found in the literature (Herron and Langway, Reference Herron and Langway1980; Fujita and others, Reference Fujita, Matsuoka, Ishida, Matsuoka and Mae2000; Cuffey and Paterson, Reference Cuffey and Paterson2010; Hörhold and others, Reference Hörhold, Kipfstuhl, Wilhelms, Freitag and Frenzel2011).
Antenna response in media and sensitivity to RIP
Antennas emit/receive electromagnetic radiation, converting the emitted/captured energy from/to electrical currents. The frequency bandwidth over which an antenna responds depends on antenna construction and geometry. The peak response (‘resonance’) is typically obtained at an in-air frequency 
$f=c/L$, with 
$L$ the characteristic length scale of the antenna itself, and should occur when the antenna currents form standing waves along the length of the antenna, implying resonant frequencies in a vacuum 
$f_{res}=mc_0/2L$, with 
$m$ some integer. For a center-fed fat dipole antenna (such as those used by RNO-G), the resonant current amplitude must be maximal at the feed-point connecting the antenna to the cable conveying the signal, corresponding to 
$m$ odd. The antenna used in this work has two such in-air resonances in the sensitive frequency band, at 
$\sim 250\,\mathrm{MHz}$ (
$m=1$) and also 
$\sim{750}\,{\mathrm{MHz}}$ (
$m=3$).
The inductive, capacitive and resistive characteristics of antennas determine their complex impedance 
$Z_A$ at a given frequency. The mismatch between the impedance of an antenna and the impedance of the input/output at the feed-point (where the antenna attaches to the cable [‘transmission line’] conveying the signal in/out of the antenna to/from a measurement device such as a Vector Network Analyzer [VNA]) determines how well it will radiate/respond at a given frequency. The complex ‘reflection coefficient’ 
$\Gamma$ is the ratio of the amplitude of the reflected wave and the incident wave in an antenna (similar to the general case of electromagnetic radiation being reflected/transmitted at a dielectric boundary) and provides a measure of the frequency-dependent antenna response; 
$\Gamma$ is defined as
\begin{equation}
\Gamma = \frac{Z_A - Z_0}{Z_A+Z_0},
\end{equation}with 
$Z_0$ the transmission line impedance (usually 
$50\Omega$ for standard coaxial cable and purely real). The magnitude of the complex 
$|S_{11}|$ parameter is related to the magnitude of the reflection coefficient via:
\begin{equation}
|S_{11}|(dB) = 10\log_{10}(|\Gamma|),
\end{equation}so 
$|S_{11}|$=0 dB corresponds to complete power reflection and therefore zero transmission; conversely, 
$|S_{11}|\to-\infty$ dB corresponds to perfect radiation at a given frequency (For the remainder of this article, 
$S_{11}$ denotes the magnitude of the 
$|S_{11}|$ parameter and not the complex value).
The medium into which an antenna is embedded can be characterized by its relative permittivity 
$\varepsilon_r$, relative permeability 
$\mu_r$ and conductivity 
$\sigma$. The in situ media complex impedance 
$Z$ then varies as:
\begin{equation}
Z = \sqrt{\frac{\mu}{\varepsilon}} = \sqrt{\frac{\mu_o \mu_r}{\varepsilon_0\varepsilon_r}} = Z_v \sqrt{\frac{1}{\varepsilon_r}} = \frac{Z_v}{n}
\end{equation} where 
$\mu_0$ is the relative permeability of the vacuum, 
$\varepsilon_0$ is the relative permittivity of vacuum and 
$Z_v$ is the free-space vacuum impedance (377
$\Omega$). We assume that the media of interest to this work are nonmagnetic (
$\mu_r = 1$). Qualitatively, this permittivity dependence implies that an antenna with measured in-air frequency-dependent impedance 
$Z(\omega)$ evolves to an in-medium impedance as: 
$Z(\omega)\to (Z/n)(\omega/n)$, separately for both the real and imaginary components of 
$Z$. In our case, for which a vertically oriented antenna is deployed into a dry borehole, the antenna is not completely embedded in a single medium. Instead, it is immersed in a medium with an effective index of refraction 
$\textit{n}_{\text{eff}}(z)$, which is a weighted combination of the 
$n\approx$1 air in the borehole and 
$\textit{n}_{\text{firn}}$ of the surrounding ice, so 
$\textit{n}_{\text{air}}\leq \textit{n}_{\text{eff}}\leq \textit{n}_{\text{firn}}$. Additional secondary effects like temperature, moisture content and impurities are assumed to be negligible for this measurement.
In this paper, we conduct experiments to demonstrate that measurements of this resonant frequency shift can be used to infer the local RIP.
Summit station measurements
Initial measurement campaign (2024)
A 
${350}\,\mathrm{m}$ deep, 98
$\pm$1 mm diameter ice hole drilled at Summit Station by the Ice Drilling Program (IDP, based at UW, Madison) using the ‘700 drill’ in June/July 2024 for PI N. Saltzman. The location of the IDP-drilled hole, relative to other landmarks at Summit Station, is indicated in Fig. 2.
Location of the 2024 IDP hole at Summit Station.
This borehole provided an opportunity to collect antenna impedance data as a function of depth. For measuring the antenna response to a depth of 
$100\,\mathrm{m}$, we deployed a fat vertically polarized antenna (VPol; Fig. 3)
(a) KU-VPol antenna schematic showing dimensions of antenna employed for data-taking and feed-point detail. An N-connectorized cable is threaded through the left cylindrical chamber to the central feed-point and then connected to the VNA. This antenna is similar in design to VPol antennas used for main RNO-G data-taking. (b) CAD model of the same dipole antenna, illustrating the central collar used to maintain the desired spacing between cylindrical halves, and also showing the endcap axial stabilizer fins, used for 2025 measurements.
constructed at the University of Kansas (KU), based on the previous design used for transmitter measurements taken at the South Pole Ice Core Experiment borehole in December 2018 (Allison and others, Reference Allison2020). Typical of dipole antennas, the in-air antenna beam pattern favors reception in the mid-plane (aka ‘boresight’, and transverse to the long axis of the antenna) with a maximum gain of +3 dBi, falling off approximately as 
$\cos\theta$ with elevation angle 
$\theta$, and approximately uniform in azimuth (
$\phi$), as illustrated in Fig. 4.
Beam pattern of dipole antenna used for data-taking, showing elevation (a) and azimuthal (b) gain. Measurements were taken in the KU anechoic chamber.
As the antenna was lowered, reflection coefficient data were collected at 
$1\,\mathrm{m}$ increments using a handheld FieldFox N9913B VNA. For these measurements, the VNA frequency resolution was set to 1 MHz. Although the highest possible frequency resolution (10,000 points over a 1 GHz frequency range, corresponding to 100 kHz resolution) was desirable, the limited battery lifetime of the N9913B FieldFox VNA imposed a practical limitation on frequency resolution, since data-taking time scales with the number of frequency bins. We note that the intrinsic frequency resolution of this instrument is quoted in the manual to be less than 1 Hz (https://www.keysight.com/us/en/assets/7018-06516/data-sheets/5992-3702.pdf). After calibration (done prior to every set of measurements), the FieldFox N9913B VNA has a systematic frequency uncertainty of 
$\pm$0.5 ppm (parts per million), such that, for example, at 1 GHz, the actual frequency uncertainty is 500 Hz. Within the limits of our resolution, we verified that the VNA reproduced the same 
$S_{11}$(frequency) profile in 16 separate S11 measurements, with no measurable scatter.
The manufacturer’s claimed depth accuracy of the GV530 winch (www.geovista.co.uk) used to lower the antenna in 2024, for which the antenna signal cable was snaked around the winch cable, is 
$1\,\mathrm{cm}$ (with no hysteresis), so error and uncertainty in the height of those measurements are considered to be negligible.
Figure 5 displays some of the recorded 2024 
$S_{11}$ profiles. We note an apparently anomalous response in the 2024 
$S_{11}$ data appearing at 
$79\,\mathrm{m}$ and then gradually decreasing over the next 
$5\,\mathrm{m}$ (discussed below).
Measured 2024 
$S_{11}$ at various depths in ice, illustrating an anomalous response observed at depths near 80 m (orange). Note the abrupt change in the shape of the resonance.
Additionally, at depths shallower than 
$20\,\mathrm{m}$, we observe significant scatter in the 
$S_{11}$ data points, indicating either some instrumental effect or true density fluctuations in this shallow depth interval. We also observe coherent reflections from small impedance mismatches, becoming increasingly visible at high frequencies as a low-amplitude ripple in the 
$S_{11}$ data, and particularly prominent at high frequencies. The period of this ripple was observed to scale, as expected, on the length of cable used for the measurement.
Improved (2025) measurements
The 2024 measurements identified one clear systematic uncertainty (possible tilt of antenna in the ice hole) and (at least) two unresolved science questions, namely the source of the scatter in the depth interval near 80 m and also shallower than 20 m. To investigate these questions further and also extract an RIP from the higher-frequency 
$m$=3 resonance, a second measurement campaign was conducted at Summit Station in May 2025.
Figure 6 shows the antenna at the outset of the 2025 measurements; evident in the picture is the large clearance of the borehole (approximately twice the nominal diameter) at the throat (unfortunately, a mechanical caliper was not available for mapping the throat profile). In 2025, the antenna was lowered by hand to eliminate any possible RF emissions due to a proximal generator or winch, using lower-loss coaxial cable (to improve the Signal-to-Noise Ratio at higher frequencies), to a depth of 
$85.5\,\mathrm{m}$ with an estimated uncertainty of 
$10\,\mathrm{cm}$. Spacers introduced in 2025 to the endcaps of the antenna attempted to address the dominant uncertainty in the 2024 data analysis arising from unknown axial alignment (discussed in more detail below). VNA sweep averaging (
$N_{avg}$=16 sweeps per data record) in 2025 also improved the visibility of the higher-frequency 
$m$=3 (
$f_2$) harmonic, in addition to the 
$m$=1 
$c_0/2L$ (
$f_1$) fundamental. For the 2025 dataset, measurements were taken as the antenna was both raised and lowered, providing cross-checking of the two datasets. Data-taking parameters for the 2 years are compared in Table 1.
Pictures of the dipole antenna (2025) at the IDP borehole. (a) In-air 
$S_{11}$ measurements, with Summit Station in the background. (b) Shows antenna at outset of data-taking; (c) illustrates large clearance of borehole at the surface, relative to antenna. Also shown is the wooden brace used (2025 only) to ensure vertical orientation of antenna and also for strain relief (assumed to have no effect on measured 
$S_{11}$ after calibration). For 2024 measurements, the signal cable was co-lowered with the strength member cable from the GeoVista GV530 winch (not shown), which was mounted atop a stationary wooden platform for stability.
Data-taking parameters for the 2024 vs 2025 datasets.
Extraction of 
$S_{11}$ resonant frequencies from raw
$S_{11}$ data
We extract, from data, the resonant frequency of the antenna in the IDP hole as a function of depth 
$f_r(z)$, by fitting the 
$S_{11}$ distribution at each depth to either a double-Gaussian (2024 data), defined as
\begin{equation}
P(f) = A_1\exp\bigg(\frac{-(f-\eta_1)^2}{2\sigma_1^2}\bigg)+A_2\exp\bigg(\frac{-(f-\eta_2)^2}{2\sigma_2^2}\bigg)
\end{equation}with 
$A_{1,2}$ the amplitude of the Gaussian, 
$\eta_{1,2}$ the mean and 
$\sigma_{1,2}$ the standard deviation, or a double-Lorentzian (2025 data):
\begin{equation}
P(f) = \frac{A_1}{(f-\eta_1)^2+\gamma_1^2}+\frac{A_2}{(f-\eta_2)^2+\gamma_2^2}
\end{equation}with 
$\gamma_{1,2}$ is the width.
For 
$m$=1 (2024 data), the double Gaussian was fit to 
$S_{11}$ over the range 1–
$300\,\mathrm{MHz}$. Sample fits are shown in Fig. 7.
Double Lorentzian fit to 2025 
$f_1$ 
$S_{11}$ profile (left) and double Lorentzian fit to 2025 
$f_2$ 
$S_{11}$ distribution using data taken at a depth of 
$13\,\mathrm{m}$ (selected arbitrarily).
Additional 
$f_{res}$ estimators were examined (tracking the 
$S_{11}$ minimum value, calculating the minimum of a weighted sum, or using an asymmetric Laplacian form, e.g.) and found to yield similar results, though with less fitting stability or poorer convergence characteristics compared to the double Gaussian/Lorentzian. The extracted resonance values 
$f_{res}(z)$ for our datasets are shown in Figure 8.
Raw extracted resonant frequencies, as returned from double-Gaussian/double-Lorentzian fit, for each indicated dataset, plotted against antenna depth.
In the 2025 data, there is a clear 2–3 MHz systematic reduction of the resonant values obtained during the (later) upwards-going (designated 
$f_{\uparrow}$) compared to the previous downwards-going data 
$(f_{\downarrow})$ for both 
$f_1$ and 
$f_2$, particularly for 
$z \lt 50\,\mathrm{m}$. Several possible causes for this discrepancy were considered. In the field, accumulated snow was observed both above and also between each of the endcap stabilizing fins when the antenna was raised from the borehole after data-taking, presumably due to the antenna scraping the side of the borehole while being lowered/raised. Such snow could affect the measured resonance by raising the effective RI around the antenna (by displacing air with snow), resulting in a lowered measured resonant frequency. Additionally, there could be VNA temperature dependence (expected to become more obvious with time), as the 2025 data were taken while the ambient air temperature was 
$-25^{\circ}$C, with a wind chill of 
$-60^{circ}$C. Both of these could bias later VNA readings towards lower values (although we were unable to intentionally reproduce a systematic, temperature-dependent droop of the VNA response in the lab). We empirically account for this by taking the mode of a distribution of 
$(f_\uparrow(z)-f_\downarrow(z))$ (Fig. 9) and applying that offset to the 
$f_\uparrow$ data.
Difference of 
$(f_{\uparrow})$ and 
$(f_{\downarrow})$ for 
$f_1$ (a) and 
$f_2$ (b) respectively (2025 data-taking). The mode bin is labeled.
Analysis
Measurement strategy
From our reflection coefficient measurements, we calculate the variation of the 
$m$=1 fundamental resonant frequency (denoted 
$f_1$) with depth for data taken in both 2024 and 2025, as well as the 
$m$=3 second resonance (first harmonic, and denoted 
$f_2$) for 2025, to derive the desired RIP. We must therefore translate our 
$S_{11}(z)$ measurements (made in a medium with local effective RI 
$\textit{n}_{\text{eff}}(z)$) into an ice 
$n(z)$ profile; this requires that our frequency dependence be re-formulated as an RI dependence, that is, 
$S_{11}(z)\to S_{11}(n(z))$, which we can then invert to obtain the desired 
$n(z)$. We expect the fundamental resonance (
$m$=1) of the antenna to scale as 
$f_{m=1} \sim c_0/(2\lambda_{m=1})$, where 
$\lambda_{m=1}$ is the 
$m$=1 fundamental wavelength. Since 
$\lambda_1$ is set by our (fixed) antenna length, and since 
$c$ scales as 1/
$n$, the resonant frequency should also scale as 1/
$n$. We correspondingly assume the functional form 
$f(n)=a/(b+n)$, with the constant 
$a$ primarily determined by the antenna length, and the constant 
$b$ allowing for the possibility of an offset independent of 
$n$.
Determination of RIP from frequency data
To extract the constants 
$a$ and 
$b$ in our parameterization for 
$f_{res}(n)$, we plot the measured resonant frequencies against the RI implied by the independent density data, described in Windischhofer (Reference Windischhofer2024). Four points (two at the endpoints of the interval, and two intermediate points) are selected over the depth interval (22–80 m) for which the frequency data are ‘smooth’; we find that our 
$f(n)=a/(b+n)$ ansatz gives a reasonable match to the observed NPM density data (Fig. 10). Having calibrated our 
$n(z)$ functional form, we can then extrapolate into the shallow firn and quantify RI fluctuations, independent of the calibration dataset; we later check against gravimetric measurements made on cores taken from the same hole.
Fits to 
$f(n)$ for 2024 
$f_1$, 2025 
$f_1$ and 2025 
$f_2$ datasets (detailed below), using four depth points selected from the ‘smooth’ depth region.
Table 2 summarizes the values of 
$a$ and 
$b$ we extract from our fits to 
$f(n)=a/(b+n)$. To assess the goodness-of-fit of our 
$f(n)\sim 1/n$ ansatz, the second column from the right tabulates the average residual of our fit (‘
$ \lt \delta \gt $’), relative to our data points (i.e., 
$\sqrt{\Sigma_i(f_{measured,i}-f_{fit})^2}$/N), where N is the number of points used in the fit, 
$f_{measured,i}$ is the fit resonant frequency for each depth value and 
$f_{fit}$ is the value of frequency, for a given depth value returned by the fit. The right-most column shows, for comparison, the deviation obtained using a linear fit 
$f(n)=a+bn$, which yields values typically larger by 
$\sim$20%, indicating preference for the inverse, rather than linear, scaling.
Fitted values of (a,b) and goodness-of-fit metrics. Last column verifies that a linear, rather than inverse dependence of frequency with RI results in a poorer fit.
To estimate our systematic error due to our selection of the four fit points, we also compared with a fit using all the data points and tabulated the difference in the extracted RIP values. Deviations are typically 0.01 units of 
$n$.
Using the extracted values of 
$a$ and 
$b$ for each dataset, we can now translate our tabulated values of resonant frequency as a function of depth 
$f_{res}(z)$ from Summit data to obtain the 
$n(z)$ profiles shown in Fig. 11. Recently, gravimetric density data for the IDP borehole were made available to our group; this provides a powerful cross-check of our 
$n(z)$ extraction procedure. Those data, in general, show good agreement with our 2025 measurements, as shown in Fig. 11.
Overlaid 
$n(z)$ profiles for each dataset, as well as the preferred RNO-G n(z) parameterization from (Windischhofer, Reference Windischhofer2024). ‘IDP Borehole’ data has slightly different binning referring to density data taken from the IDP hole, converted to RI assuming the linear scaling 
$n(z)=1+0.8485\rho(z)$.
Restricted to the shallow depth region (
$z \lt $20 m), we observe that the density-derived IDP hole RIP visually tracks our 
$S_{11}$-derived RIP (Fig. 12a), averaging 
$f_{1,\uparrow}$, 
$f_{1,\downarrow}$, 
$f_{2,\uparrow}$ and 
$f_{2,\downarrow}$, indicating that both techniques are sensitive to real density fluctuations in the shallow firn.
(a) Comparison of 
$n(z)$ profile based on IDP-borehole density data with 
$n(z)$ profile based on 
$S_{11}$ IDP-borehole measurements, for shallow depths. (b) Deviation between the indicated 
$n(z)$ profile, relative to the RNO-G default parameterization. (c) One-dimensional projection of residuals in the previous plot over the interval from 20–80 m; means and widths of distributions are: 
$ \lt \delta f_1(2024) \gt =(3.5\pm1.4)\times 10^{-3}$, 
$ \lt \delta f_1(2024) \gt =(4.1\pm3.6)\times 10^{-3}$ and 
$ \lt \delta f_2(2025) \gt =(3.2\pm4.4)\times 10^{-3}$. We interpret the general consistency with zero offset as lending support to the assumption of inverse RIP scaling with frequency, used to extract the 
$n(z)$ profile.
We have performed additional cross-checks of our functional parameterization for 
$f(n)$. Panel (b) of Fig. 12 shows the deviation between our extracted 
$n(z)$ values relative to the RNO-G parameterization; panel (c) displays the distribution of deviations for the three datasets considered. Overall, in the ‘good’ region (20–80 m), we observe that the deviations, relative to the 
$n(z)$ parameterization, are flat (and therefore unbiased in depth); the magnitude of the mean of the deviations is of order 0.003 units of 
$n$.
For comparison, we have tabulated the NPM density data taken in the vicinity of Summit and translated to RI for comparison (Fig. 13). The 2024 IDP density data are also overlaid. The profile shapes generally track each other; we note significant variations in the calculated RIP over sub-km spatial scales, comparable to the separation between the 350 m IDP hole and the nominal Summit Station NPM data-taking site, underscoring the importance of local measurements. Similar to the trend observed in our 2025 
$n(z)$ profile, the IDP borehole density data overshoots the NPM data in the shallow depth regime.
Comparison of 
$n(z)$ profiles implied by Hawley& Morris (private communication) and (Morris and Cooper, Reference Morris and Cooper2003; Hawley and Morris, Reference Hawley and Morris2006) 
$\rho(z)$ NPM data taken in the vicinity of Summit Station (retaining original site-naming convention; these data are one input to the 
$n(z)$ parameterization favored by RNO-G). Also included are density data taken at the IDP hole in 2024. Distances indicated in the key are relative to the Summit Station main base. At shallow depths, we note significant scatter between sites separated by less than 1 km laterally.
Consistency between the 2024 and 2025 datasets and uncertainties
We have attempted to quantify the level of agreement between the August 2024 and May 2025 data samples. Although they are consistent in the intermediate depth interval (20–80 m), they conspicuously deviate at depths of 80–85 m, for which the 2024 
$f_1$ data ‘overshoot’ the density-derived RI, and in the shallow firn (depths less than 22 m). Whereas the 2024 data-taking employed a (large metal) winch at the top of the hole, with signal and winch cable both in the hole, the 2025 data-taking used a single coaxial cable, eschewing use of the winch or ancillary cables. Perhaps most importantly, there was no attempt to axially center the antenna in the initial 2024 data-taking.
The comparison is also complicated by the 2–3 MHz systematic offsets between different data samples, as well as the possible effects of additional 1 year of snow accumulation at the IDP hole site (averaging 24 cm water-equivalent yearly, at Summit Station) or densification effects resulting from human activity, including heavy machinery, at the drill site in the intervening time period). Ignoring any such time-dependent effects, a direct Pearson’s coefficient calculation gives a correlation coefficient of 0.537. Alternatively, we can allow for a depth offset (but less than 1 m) between the two data samples by calculating the binary coherence between the shapes of the two 
$n(z)$ curves, as measured by whether the 
$(i+1)$th data point is ‘coherently’ higher/lower than the 
$i$th data point. We find 17 cases for which a point-to-point increase/decrease in the 2024 data coincides with a correlated increase/decrease in the 2025 data, compared to three instances where they are anticorrelated, giving a binary probability in excess of 
$ \gt $99% correlation in the shallow region. Importantly, we find that the 
$S_{11}$ fluctuations observed in the shallow depth interval are reproduced for the 2025 up-going vs down-going datasets, indicating that these are real density variations and not instrumental effects.
To investigate the effect of the 2025 stabilizer, we took dedicated runs in May 2025, with and without the stabilizer, down to 20 m, observing a narrower distribution with the stabilizer (Fig. 14) in place. We note that the diameter of boreholes is typically largest near the surface – we measured a throat diameter of 21.8 cm of the borehole at the surface in 2025, indicating that off-axis effects, particularly near the surface, could be severe. Moreover, simulations indicate that, if the antenna were inclined on the hole wall, we should observe a systematic shift in resonant frequency (and therefore RI).
Comparison between 
$f_{\uparrow}-f_{\downarrow}$ for 
$f_1$ (a) and 
$f_2$ (b) frequencies extracted for runs (depths shallower than 20 m) using antenna spacers for centering/axial stability (2025 data configuration; green) vs runs without stabilizers (2024 data configuration; orange). We observe increased consistency between the up/down data samples with spacers in place.
For 2024, each depth only had one measurement, rendering it unfeasible to estimate a true statistical uncertainty. Instead, we estimate the systematic error from the measuring device by scaling the uniform errors of the 
$S_{11}$ profile until the fit has a reduced chi-squared 
$\chi^2/\text{ndof} = 1$ (A fuller description of the chi-squared rescaling (standard in particle physics) can be found in the literature (Tanabashi and others, Reference Tanabashi2018; Cowan, Reference Cowan2019, Introduction (pp. 1–19))). The Hessian matrix error on the returned fit 
$\eta_1$ value at a given depth was used as the VNA systematic uncertainty. This is only an approximate uncertainty as the double Gaussian parameterization is empirical, and not a first-principles description of the underlying 
$S_{11}$ trace distribution. Therefore, the uncertainties extracted from the fits should only be viewed as a proxy for the true statistical uncertainties of taking multiple measurements at the same depth. The uncertainty in the parameterization used to determine the depth-dependent resonant frequency can be assessed by considering the variation between the extracted 
$n(z)$ profile using three fitting parameterizations: i) the double Gaussian or ii) Lorentzian described earlier (our default) and iii) the Doniach–Sunjic line shape (Doniach and Sunjic, Reference Doniach and Sunjic1970), which is found to empirically give a reasonable match to our 
$S_{11}(f)$ line shapes. The deviation between the three fitting schemes has a center inter-quartile range of 
$\Delta n=0.0064$, indicating reasonable agreement. For the 2025 dataset, the magnitude of errors is estimated directly from the discrepancy between up/down, at each data point.
There is a systematic uncertainty from the measurement device itself due to limited instrumental resolution. Other systematic errors inherent in the network analyzer include possible calibration drift introduced by any bends in the cable and or any temperature variation during data-taking (the VNA is nominally rated for operation between 
$-10^{\circ}$C and 
$+55^{\circ}$C).
We tabulate the considered systematic uncertainties in Table 3.
Systematic uncertainties and the estimated effect on the measured resonant frequency.
Conclusions, discussions and future work
Understanding and quantifying ice properties such as the RI are critical inputs to polar UHEN experiments. Three different accuracy ‘scales’ are relevant for neutrino measurements—a) neutrino ‘identification’ can be done by ignoring the RI variation with depth and using constant 
$n(z)$, since the primary parameter that separates radio signals produced by neutrino interactions from backgrounds is a simple up/down criterion—the former are upcoming, while backgrounds (primarily anthropogenic, but also including ‘natural’ backgrounds such as those generated by the tribo-electric effect or down-coming cosmic rays) arrive from the surface, b) once identified, translating a statistically large sample of neutrinos into a flux requires an estimate of the neutrino ‘effective volume’, defined as the equivalent ice volume over which neutrino interactions can be measured. For low detectable neutrino energies (10–100 PeV), the extent of the shadow zone dominates the estimated effective volume (
$V_{\text{eff}}$), with variations of 1% in RIP translating into uncertainties in 
$V_{\text{eff}}$ an order-of-magnitude larger (Ali and others, Reference Ali2024). Uncertainties in the 
$n(z)$ profile therefore directly translate into weaker upper limits on the neutrino flux at this energy scale. At ‘high’ neutrino energies, 
$V_{\text{eff}}$ is determined by ice attenuation, so RI systematic uncertainties are less important. For neutrino ‘astronomy’, for which signal reconstruction depends on precise ray tracing, 1% RI uncertainties again translate into 10
$\times$ larger uncertainties in arrival direction.
In this paper, we have attempted to correlate antenna response with the index of refraction profile as a function of depth. We find agreement between the shallow firn fluctuations implied by our IDP borehole 
$n(z)$ profile with direct density measurements, taken from the same borehole. Compared to other techniques that can be used to infer an RIP, NPM (e.g.) is statistically superior to what can be achieved with antenna resonance. The advantage of the approach outlined here is speed (the 2025 data-taking required approximately 30 min of run time), does not require a winch or dedicated power source, uses hardware readily accessible to the UHEN experiments (antennas and a VNA), can be easily integrated into the RNO-G field plan and should be sensitive to RIP differences between two sites laterally separated by the 1–2 km typical spacing between UHEN receiver stations.
Future measurements
Although UHEN experiments favor broadband antennas to maximize sensitivity to the impulsive, ns-scale signals produced by neutrinos interacting in-ice, RI measurements using this technique favor ‘ideal’ (long, thin wire) narrow-bandwidth dipole antennas with high Q-values deployed into similarly narrow boreholes. The Stampfli drill operated by the IDP, for example, produces 57 mm diameter boreholes to a depth of 20 m in 3 hours. Such a measurement would thereby minimize systematics due to air content in the hole as well as fitting uncertainties inherent in fitting 
$S_{11}(\omega)$ profiles. Even in such an ideal case, however, an irreducible systematic error due to the fact that one must sample over meter-scale wavelengths, and therefore average over smaller-scale fluctuations, remains. In 2025 and 2026, the RNO-G experiment will systematically obtain additional 
$S_{11}$ from 
$\approx$10 newly drilled holes at the Summit Station site, offering a considerably extended data sample for similar studies in the future.
Acknowledgements
We thank the support staff at Summit Station for making RNO-G possible, as well as Jay Johnson and the UW-based IDP team for drilling the icehole from which the primary measurements in this paper were derived. We are deeply indebted to the KU Physics and Astronomy Machine Shop, and particularly Scott Voigt and Mark Stockham, who designed and constructed the custom antenna used for recording the primary data essential to this measurement, as well as the spacers employed for antenna stabilization in 2025. We also acknowledge our colleagues from the British Antarctic Survey for building and operating the BigRAID drill for our project. We are indebted to Murat Aydin, Eric Saltzman and John Patterson for supplying the density vs. depth measurements for the 350-m borehole, which provided an essential reference standard. Jay Johnson provided essential details on the specifications of the ice drill, as well as the borehole, drilled by IDP during the summer of 2024. We also thanks Liz Morris and Bob Hawley for insightful conversations, and for availing us of their field experience, and also the NPM data to which we compared our extracted n(z) profile, as presented in this manuscript.
We would like to acknowledge our home institutions and funding agencies for supporting the RNO-G work; in particular the Belgian Funds for Scientific Research (FRS-FNRS and FWO) and the FWO programme for International Research Infrastructure (IRI), the National Science Foundation (NSF Award IDs 2112352, 2111232, 2111410, and collaborative awards 2310122 through 2310129), and the IceCube EPSCoR Initiative (Award ID 2019597), the Helmholtz Association, the Swedish Research Council (VR, Grant 2021-05449 and 2021-00158), the University of Chicago Research Computing Center, and the European Union under the European Unions Horizon 2020 research and innovation programme (grant agreements No 805486), as well as ERC (Pro-RNO-G No 101115122 and NuRadioOpt No 101116890).
Open access
