Electric field reconstruction with three polarizations for the radio detection of ultra-high energy particles

We used the ZHAireS simulation package [19] to simulate the radio emissions induced by EAS. We created a Monte Carlo (MC) dataset covering zenith angles $\theta_{0}$ from 63.0^∘ to 87.1^∘ in equidistant steps of $\log_{10}(1/\cos(\theta_{0}))$, focusing specifically on inclined EAS. The azimuth angle $\phi_{0}$ spans from 0 to 180^∘, with increments of 45^∘. The primary particle energy $E_{\mathrm{shower}}$ ranges from 0.126 EeV to 3.98 EeV in logarithmic steps of 0.1. The simulation site is located in Dunhuang, a radio-quiet area in China. The geomagnetic field at this site has a strength of 56$\mu$T and an inclination angle of 61 ^∘. For the atmosphere, we use the extended Linsley standard atmosphere, in which an exponential index of refraction model was used, with scale height of 8.2km and sea level refraction index of 1.000325. Our simulation set is evenly divided between proton and iron primaries. A total of 4160 EAS events were generated for this analysis.

We employed a star-shaped configuration on the ground comprising 160 antennas in the simulation, organized into 8 arms with 20 antennas each. The ring spacing was adjusted to align with the Cherenkov angle, ensuring that the peak amplitude was located centrally within the entire simulated footprint. This adjustment aims to ensure robust footprint coverage, even under aggressive trigger conditions.

The components of the simulated electric field correspond to three polarization directions: $\mathrm{E_{x}}$ (South-North), $\mathrm{E_{y}}$ (East-West), and $\mathrm{E_{z}}$ (Down-Up). We selected a time bin of 0.5 ns within a 1000 ns time window, which provides a frequency resolution of 1 MHz after performing a Fourier transformation.

When using only horizontally aligned antennas, a significant portion of events coming close to the horizon cannot be detected, as a significant part of the signal will be vertically polarized. To achieve full sky coverage, it is essential to include an additional antenna arm oriented vertically [40], or have two antennas that are partly sensitive to the vertical component as the case for the SALLAs in Auger-RD [41].

This study tests two types of antennas to reconstruct the electric field with three polarizations: a simple 3-arm dipole antenna, which has symmetrical oscillators for the horizontal arms, and a more complex 3D bow-tie-shaped HORIZON antenna [22, 42]. The dipole antenna serves as a baseline for developing the reconstruction method, while the HORIZON antenna is used for verification and performance evaluation for a realistic setup. Antenna responses were simulated using Ansys HFSS [43], with all simulations incorporating ground effects to ensure accurate representation of real-world conditions. The antenna response patterns are shown, for a representative frequency, in Appendix A.

The open-circuit voltage is a convolution of equivalent length (antenna response) as shown in Figure 1 and 2, and electric field:

where $V(t)$ represents the measured voltage, $\vec{H}(t)$ is the vector equivalent length (VEL), $\vec{E}(t)$ is the electric field, and $*$ denotes the convolution operation. For simplification purposes, calculations are often performed in the frequency domain through Fourier transformation, yielding the expression [47]:

where $\mathcal{V}(f)$, $\vec{\mathcal{H}}(f)$ and $\vec{\mathcal{E}}(f)$ are voltage, VEL, and electric field respectively in the frequency domain after Fourier transformation.

Reconstructing radio signals is challenging due to the complexity of various background noise sources present at the site. Most of the radio emission power is concentrated in the frequency range below 200 MHz. To mitigate the impact of human activities such as radio telecommunications and navigation, radio detection experiment sites are chosen to be located far from urban areas. While thunderstorms, solar flares, subterranean pre-seismic movements, and other intense radio activities can generate complicated background noise, they are only active at certain times and are not included in this work. In this scenario, the galactic background is considered to be the predominant noise source in the specified frequency band [48]. This type of noise, which marked the birth of radio astronomy, was first identified by Karl Jansky and originates from the diffuse synchrotron radiation within the Galaxy [49].

To model the background noise spectrum in the range of [30, 200] MHz, the LFmap model [50] is utilized. This model includes various components such as the Cosmic Microwave Background (CMB), isotropic elements, Galactic sources, and extragalactic sources. The sky temperature follows a power law in this context:

The brightness can be derived from the sky temperature using the Rayleigh-Jeans approximation: $B_{\mathrm{f}}=\frac{2k_{\mathrm{B}}f^{2}T_{\mathrm{sky}}}{c^{2}}$, where $k_{B}$ is the Boltzmann constant, and $c$ is the speed of light. The power of sky noise is then expressed as:

where $t^{\prime}$ represents the local sidereal time (LST), $A_{e}(\theta,\varphi,f)$ is effective receiving area of the antenna (obtained from the antenna gain), and the antenna efficiency $\eta_{A}=1$. The effective area is given by $A_{e}=\frac{\lambda^{2}}{4\pi}G$, where $G$ is the antenna gain, G=1 for a lossless resonant antenna of impedance matching the load. The left panel of Figure 3 shows this background power spectrum in East-West polarization.

Here, we do not consider any specific antenna design but assume a typical antenna impedance of $Z_{0}=50\Omega$, corresponding to an ideal impedance match. In practice, antenna impedance varies with frequency and may deviate from this idealized assumption. However, this simplification provides a consistent noise level, ensuring a reliable baseline for comparison despite potential discrepancies with real-world scenarios. With impedance matching, half of the power is transferred to the receiver. Consequently, the voltage received from the background can be expressed as:

The right panel of Figure 3 shows the RMS of the Galactic noise in the time domain as a function of the LST in East-West polarization.

Since the voltage spectrum of the noise is known and can be approximated as white noise, we generated a uniformly distributed phase that was added to this spectrum. By performing an inverse Fast Fourier Transform (iFFT), we obtained the noise trace in the time domain, which was then added to the signal’s voltage trace.

Thermal noise in the electronics of the Radio Frequency (RF) chain also contributes to the overall noise. However, this noise is added independently at each polarization after the open circuit voltage. Since our focus is on the reconstruction from open circuit voltage to electric field, electronic noise is beyond the scope of this article.

When the SNR is large enough in all the traces in the full frequency band, and given the signal’s arrival direction, two orthogonal horizontal polarizations are sufficient to give a good reconstruction of the electric field. This has been implemented in most radio experiments [47, 51]. The electric field can be reconstructed by inverting Equation 3.2 using any two voltage polarizations, expressed as $\vec{\mathcal{E}}=\vec{\mathcal{H}}^{-1}\mathcal{V}$. In this approach, noise is considered to be small enough to assume that $\mathcal{V}_{\mathrm{oc}}(f)$ is equal to $\mathcal{V}(f)$.

In the absence of noise, this method yields an exact solution for the electric field. With the inclusion of an additional polarization, the electric field can be reconstructed in three different combinations, as demonstrated first in LOPES 3D. Since different polarizations are sensitive to varying incident angles, LOPES 3D employed a strategy of weighting the reconstructed electric field based on the zenith angle [44].

Although the $\hat{r}$ component is small, it is non-zero, making it possible to solve the inverse matrix in Equation 3.1. This allows us to include the third polarization, which is expected to improve the signal by providing additional information.

We compare this approach with the conventional 2 polarizations method using a simple dipole model in the 30-80 MHz frequency band. The example we present shows voltage traces with an SNR of approximately 6 in each polarization. Here, the SNR is defined as $S_{\mathrm{peak}}$/RMS_noise, where signal S is the maximum value of the Hilbert envelope of the voltage, and the root mean square (RMS) of noise is associated with the standard deviation in a 250 ns noise window located 500 ns after the signal peak. Figure 4 illustrates the signals in the three polarizations.

The reconstruction results for this example are shown in Figure 5, and Table 1 provides a quantitative comparison in the time domain through normalized cross-correlation values between the reconstructed and simulated (Monte Carlo true) traces. We take the normalized cross-correlation from [52] which is used to check the agreement of interpolated pulse shape and is defined as:

where $f(t)$ represents the simulated electric field trace, $g(t)$ represents the reconstructed trace, and $t_{\mathrm{max}}$ is the trace length. And $\tau$ is the time mismatching between these two traces, in our case $\tau=0$.

For the $\varphi$-component of the electric field, which carries a relatively stronger signal, both the two-polarization and three-polarization reconstruction methods produce similar results. As shown in the upper right panel of Figure 5, the simulated and reconstructed traces overlap closely. In the frequency domain (lower right panel), while the three-polarization method introduces spikes, it still achieves a normalized cross-correlation of around 0.80, slightly underperforming the two-polarization method. For the $\theta$-component (left panel), which has a weaker signal, the differences between the methods become more evident. The three-polarization method shows fluctuations at specific frequencies, whereas the two-polarization method exhibits larger deviations toward the end of the spectrum. This is reflected in the negative cross-correlation value for the two-polarization approach, indicating greater errors and a noticeable peak time shift observed in the time-domain trace.

The spikes in the spectrum for the three-polarization reconstruction, as shown in the bottom panel of Figure 5, are due to the very small values of the $\hat{r}$-component, which cause divergence at certain frequencies.

To eliminate the divergence at certain frequencies in the three-polarization reconstruction, we employ a robust least squares estimation method. This approach enables a $\chi^{2}$ minimization process without requiring any prior assumptions about the electric field spectrum, thereby enhancing the reliability of the reconstruction. By accounting for noise and minimizing the influence of small $\hat{r}$-components, this method ensures more stable and accurate results across the frequency spectrum.

We build a $\chi^{2}$ based on antenna response, noise spectrum and measured voltage:

where $\mathcal{V}$ is the voltage in the 3 polarizations, and $\mathcal{H}$ is the antenna response shown in Equation 3.2. Since the antenna response depends on the direction of the EAS, we use the true direction from simulations in this study. The choice is based on the observation that the antenna response does not exhibit abrupt changes with direction, and the directional resolution provided by the plane wave reconstruction suffices for this case. A brief discussion on the impact of directional smearing is provided in section 4.2.3. The electric field is represented as $\boldsymbol{\mathcal{E}}=\left(\begin{array}[]{l}\mathcal{E}_{\theta}\\ \mathcal{E}{{}_{\varphi}}\end{array}\right)$. We assume that the error in this measurement is introduced by background noise. Therefore, the covariance matrix $\sigma_{\mathcal{V}}=\operatorname{diag}\left(\sigma_{\mathcal{V}1},\sigma_{\mathcal{V}2},\sigma_{\mathcal{V}3}\right)$ is given by the square spectrum of it.

Taking the derivative to minimize $\chi^{2}$:

We get an analytical solution of electric field [53]:

The global minimum does not assure a minimum on each antenna arm, but by doing the minimization processes on three polarizations simultaneously, a better overall accuracy is obtained. Additionally, as with all kinds of minimization procedures, we must make sure that the scanning algorithm of the parameter space and the sample size are adequate to drop out the uncertainty in the reconstruction. In this work, the analytical least-squares ($\chi^{2}$) method directly provides the reconstructed electric field on each antenna arm using just the measured voltage, antenna response and background noise spectrum, which is more straightforward than the other methods.

In experiments, degradation of the antenna array can occur over time, some antennas may be broken or damaged, consequently their responses may also change due to some reasons. In this case, this analytical least-squares ($\chi^{2}$) method is still reliable, since the $\chi^{2}$ minimization on each antenna is not influenced by the unreliable antennas.

In Figure 6, we present a comparison between our analytical least-squares ($\chi^{2}$) method and the matrix inversion approach for reconstructing the electric field for the dipole antenna. The time-domain traces exhibit smaller fluctuations in the non-signal time window, and the spectrum is free of spikes, both of which highlight the improved reconstruction achieved using the analytical least-squares ($\chi^{2}$) method.

To evaluate the performance of the analytical least-squares ($\chi^{2}$) method, we again employ Equation 3.5 for comparison with the matrix inversion method. We found the analytical least-squares ($\chi^{2}$) method effectively minimizes noise-induced distortions and ensures more stable results. The normalized cross-correlation values listed in Table 2 show that the analytical least-squares ($\chi^{2}$) solution yields values better than those of the matrix inversion method for both components, confirming better performance in reconstructing the electric field.

To confirm that our method is not highly dependent on the specific antenna response and frequency bandwidth, we applied it to the HORIZON antenna (see Sec. 2.2.2), which operates over a broader frequency band. Using the same simulation parameters as those for the dipole antenna example, we examined the differences in reconstruction under the same near threshold conditions. This particular signal has larger horizontal component than vertical ones, as shown in Figure 7. This comparison highlights the robustness of our method across different antenna designs and signal conditions.

The reconstruction results are shown in Figure 8. The $\varphi$-component (right panel) is well reconstructed with the analytical least-squares ($\chi^{2}$) method, while the $\theta$-component (left panel), which has weaker signals, shows significant inaccuracies with the matrix inversion method (shown in coral). The bottom panel of Figure 8 compares the spectra, providing a detailed assessment of reconstruction quality.

It is important to note that with an azimuth angle of $45^{\circ}$, the response is symmetric for horizontal polarizations, and the vertical component’s response is only sensitive to the $\theta$-component. For $\varphi$-component, there is a bump that occurs around 120MHz due to the low antenna gain at this frequency,as illustrated by the solid line for the horizontal component in Figure 2. The time-domain trace is not significantly affected, as shown by the normalized cross-correlation value being close to 1 . For the $\theta$-component, the matrix inversion method performs much worse than the least-squares approach, a large deviation appears around 60MHz, mainly due to the low antenna gain in the $\theta$-component at this frequency, where the gain is primarily derived from the vertical polarization. The fluctuations at the higher end of the spectrum result from both signal weakening and reduced antenna response at those frequencies. Although both methods exhibit deviations, the analytical least-squares ($\chi^{2}$) method provides better accuracy (see normalized cross-correlation value in Table 3).

Before the statistical analysis, we must select the events that can be distinguished from the background noise. The selection criteria for identifying simulated signals are informed by the characteristics of the simulated voltage and the background noise. The following criteria are intended to ensure that the detected signals are genuine:

• •

  SNR $>$ 5 for any of the 3 polarizations.

• •

  $t_{\text{signal}}-200$ ns ¡ $t_{\text{peak}}$¡  $t_{\text{signal}}$ + 200 ns : the time of the maximum amplitude of the trace, $t_{\text{peak}}$, must fall within $\pm 200$ ns of the time the signal is expected to be on the trace.

• •

  Only the innermost 16 antennas out of the 20 simulated antennas of each arm of the star-shape pattern are used. The simulation footprint covers twice the Cherenkov angle, and antennas beyond number 16 are excluded because they do not capture significant coherent signals.

• •

  $N_{\text{antenna}}\ \ \geqslant$ 5 : at least 5 antennas must meet the above criteria for the simulation event to be selected, which is a typical requirement for cosmic ray reconstruction.

For the dipole antenna, 2332 out of 4160 simulations pass these event selection criteria. In these selected events, 79% of the antennas were triggered according to the above conditions. For the HORIZON antenna, 2749 simulations passed the selection criteria, from which 71% of the antennas were triggered. This selection process does not introduce bias with respect to the initial parameters, such as zenith angle, azimuth angle, or energy.

The Hilbert envelope offers a smooth, amplitude-tracking representation of the signal, making it especially useful for handling oscillatory or noisy data. We use the peak value of the Hilbert envelope of the reconstructed electric field for comparison between true and reconstructed traces. The relative error distribution is shown in Figure 9.

The distribution shows that the reconstruction method performs robustly for both antenna types, with the relative error distributions approximately centered, indicating reliable accuracy in reconstructing the electric field. However, differences in bias and error distribution are observed between the dipole and HORIZON antennas.

For the dipole antenna, the results exhibit a relative error of 0.03, with a mean and median bias of 0.05, which suggests a slight underestimation in the reconstruction using this method. The slight asymmetry in $\psi_{68}$ on either side further reflects this bias.

HORIZON antenna shows a smaller bias, with a mean and median of 0 indicating a closer alignment with the true signal. However, the HORIZON antenna also exhibits a slightly broader error distribution, with a standard deviation (std) of 0.04, particularly in the $\theta$-component with std of 0.12. This broader spread is likely attributed to the antenna’s response, which prioritizes sensitivity to inclined air showers. The larger frequency bandwidth of the HORIZON antenna contributes to increased variability in the $\theta$-polarization, thereby accumulating greater errors. Nonetheless, these results indicate that the analytical least-squares ($\chi^{2}$) method achieves robust reconstruction performance with the HORIZON antenna.

For both antennas, the $\varphi$-component, representing the horizontal polarization, is generally reconstructed with higher precision compared to the $\theta$-component. As shown in Figure 1 and Figure 2, the $\varphi$-component benefits from stronger signals and better antenna gain in most of the frequency band. This trend is more pronounced for the HORIZON antenna, which shows higher variability in the $\theta$-component in its broader frequency band.

In conclusion, the overall relative error of reconstruction remains small for both antenna types, despite minor variations in bias and difference in relative error. This consistency underscores the robustness of the analytical least-squares ($\chi^{2}$) method for electric field reconstruction in cosmic ray studies.

The electromagnetic component of an air shower carries the majority of the energy of the primary particle E_shower, and the electromagnetic energy E_em can be estimated directly from the radiation energy [54]. For each antenna, the radio signal pulse determines the energy deposition per unit area, called energy fluence. For an antenna array, the spatial integral over the energy fluence measures the radiation energy from the primary particle in the working frequency band of the antennas. Hence, we can reconstruct the primary energy from the reconstructed electric field.

The energy fluence is defined as the integral of the squared electric field within the signal time window, which represents the total deposited energy per area accumulated at the location of the antenna.

with $\vec{E}$ as the reconstructed electric field, c the light speed in vacuum and $\epsilon_{0}$ the permittivity in vacuum.

For coherent emission, the radio pulse amplitude is proportional to the electromagnetic energy of the air shower, so $\Phi\propto|\vec{E}|^{2}\propto$ E${}^{2}_{\mathrm{shower}}$ [54, 55, 38].

Radio emission starts from the starting point of the electromagnetic interaction, carrying information at each stage of air shower development until the end of it. The typical duration of radio emission is of order several nanoseconds. We selected a 100 ns time window symmetrically centered around the signal peak and a 100 ns window at the end of the trace to estimate the fluence of the noise. The background noise fluence is subtracted from the fluence calculated for the signal window to calculate the signal energy fluence. Figure 14 displays the relative error distribution in energy fluence, which exhibits a similar trend to the Hilbert peak distribution but with more pronounced discrepancies since the error in the energy fluence accumulates from the squared electric field error. This increased deviation is attributed to the fact that noise accumulates over time during the integration process and cannot be completely filtered out.

For the dipole antenna, the median underestimation persists at 3%, which is consistent with the bias observed in the Hilbert peak analysis. While the HORIZON antenna exhibits no bias in energy fluence. Notably, the $\theta$-component exhibits worse performance in both antenna types due to the relatively weaker signal in this component.

In terms of overall precision, the standard deviation in energy fluence is approximately 6%, with the HORIZON antenna displaying a more asymmetric $\psi_{68}$ range of [-0.04, 0.01]. In contrast, the dipole antenna shows a 4% standard deviation in total energy fluence with a more symmetric $\psi_{68}$ distribution of [-0.04, 0.03]. Despite these differences, both antennas enable precise reconstruction, demonstrating the method’s effectiveness.

Figure 15 compares the reconstructed and simulated square root of energy fluence $\sqrt{\Phi}$, which is proportional to E_em. For each event, the SNRs on the antennas vary depending on the distance and viewing angle to the shower maximum. Additionally, for each antenna, the sensitivity of the 3 arms also depend on the air shower direction, leading to SNR differences in each channel. We classify events by SNR, and we found our least-squares method yields accurate energy fluence estimates at high SNRs ($>$10) for both types of antennas, while at lower SNRs, increased deviations are observed, which indicates bigger uncertainties in the reconstruction of shower energy E_shower.