Thunderstorm electric fields probed by extensive air showers through their polarized radio emission

We observe a large fraction of circular polarization in radio emission from extensive air showers recorded during thunderstorms, much higher than in the emission from air showers measured during fair-weather circumstances. We show that the circular polarization of the air showers measured during thunderstorms can be explained by the change in the direction of the transverse current as a function of altitude induced by atmospheric electric fields. Thus by using the full set of Stokes parameters for these events, we obtain a good characterization of the electric fields in thunderclouds. We also measure a large horizontal component of the electric fields in the two events that we have analysed.


I. INTRODUCTION
Lightning initiation [1] and propagation [2] are driven by the electric fields in a thunderstorm. However, performing measurements of these fields is very challenging due to the violent conditions in thunderclouds. A nonintrusive method to probe thunderstorm electric fields is through a measurement of radio emission from extensive air showers during thunderstorms [3].
When a high-energy cosmic ray strikes the Earth's atmosphere, it generates many secondary particles, a so-called extensive air shower. The dominant contribution to the radio emission from air showers during fair weather (fair-weather events) is driven by the geomagnetic field [4,5]. Electrons and positrons are deflected in opposite directions due to the Lorentz force, which results in a current perpendicular to the shower axis. As the shower develops, this current varies with altitude, thereby producing radio emission. This radiation is linearly polarized in theê v×B direction, where v is the velocity of the shower front, B is the geomagnetic field, andê denotes a unit vector. In addition, as the shower propagates, a negative charge excess builds up in the shower front due to the knock-out of electrons from air molecules by the shower particles. The variation of this charge excess gives rise to a secondary contribution to the emission [6,7]. The chargeexcess emission is also linearly polarized, but radially with respect to the shower axis. For fair-weather events, we observe a small fraction of circular polarization due to the fact that the time structures of the radio pulses emitted from the charge-excess component and those from the transverse-current component are different [8]. Since the charge-excess pulses are delayed with respect to the transverse-current pulses and they are polarized in different directions, the polarization of the total pulse rotates from one direction to the other. In our analysis, this is seen as circular polarization where the magnitude and handedness depend on the distance and the azimuth position of the observer with respect to the shower axis.
As shown in [3], due to the influences of atmospheric electric fields, intensity and linear-polarization footprints of the showers observed during thunderstorms (thunderstorm events) are different from those of fair-weather events. In this paper we show that thunderstorm events have a larger circular polarization component near the shower axis than fair-weather events. We demonstrate quantitatively that this can be explained as being due to the variation of the atmospheric electric field with altitude. The electric field changes the direction of the transverse current and thus changes the polarization direction of radio emission. The signals from the different layers are emitted in sequence when the air-shower front, progressing with essentially the light velocity, c, passes through. The emitted radio signals travel with a lower velocity than the shower front, c=n, where n is the index of refraction. Thus, near the shower axis, the pulses from the upper layers arrive with a delay with respect to the pulses from the lower layers resulting in a change of the polarization angle over the duration of the pulse, which is seen as circular polarization in the data. Therefore, the usage of the circular polarization measurements puts strong additional constraints on the structure of the atmospheric electric fields on top of the information obtained by using only the radio intensity. Since the circular polarization is due to a reorientation of the transverse current in the shower front, the circular polarization does not depend on the azimuthal orientation of the antenna with respect to the shower axis, unlike in the case for the circular polarization of fair-weather events.
In this work, we present data on circular polarization seen in the radio emission of a large number of thunderstorm events close to the shower axis as measured with the Low-Frequency Array (LOFAR) radio telescope array; see Sec. II. In Sec. III, we present a toy model to explain the cause of circular polarization of air showers measured during thunderstorms. Two reconstructed thunderstorm events are presented in Sec. IV to show that circular polarization is essential to obtain additional information about the atmospheric electric fields. Conclusions are given in Sec. V.

II. LOFAR AND DATA ANALYSIS
Data for the present analysis were recorded with the lowband antennas (LBAs) in the core of the LOFAR radio telecope [9]. Each LBA consists of two dipoles and records in the frequency range of 10-90 MHz. These antennas are grouped into circular stations. The stations are positioned with increasing density towards the center of LOFAR. The highest density is at the core where six such stations are located in a ∼320 m diameter region, the so-called "superterp." For the purpose of air-shower measurements, these antennas are equipped with ring buffers that can store the raw voltage traces sampled every 5 ns, up to 5 s. A trigger is obtained from a particle detector array, LOFAR Radboud Air-Shower Array, from air showers with a primary energy in excess of 2 × 10 16 eV [10].
The data are processed in an off-line analysis [11]. The arrival direction of the air shower is reconstructed based on the arrival times of the radio signals in all antennas. The primary energy of the air shower is estimated by using the particle detector data. The radio signal containing the pulse is received by an antenna where the signal amplitude S i is determined at 5 ns time intervals, i.e., sampled at 200 MS=s, where the sample number is denoted by i. For each antenna, the Stokes parameters, I, Q, U and V, are expressed as as derived in Ref. [12]. ε i ¼ S i þ iŜ i are the complex signal voltages, whereŜ i is sample i of the Hilbert transform of S. The summation is performed over n ¼ 5 samples, centered around the peak of the pulse. Stokes I is the intensity of the radio emission. Stokes Q and U are used to derive the linear-polarization angle and Stokes V represents the circular polarization. During the period between June 2011 and January 2015, there were 118 fair-weather events [13] and 20 thunderstorm events [3] with radio signals detected in at least four LBA stations. For comparison, the circular polarization for 20 thunderstorm events and for six fair-weather events is shown in Figs. 1 and 2, respectively. The circular polarization for fair-weather events is very small near the shower axis and increases with distance [8]. Therefore, in order to show the dependence on azimuth angle, φ, we selected those fair-weather events that have data of at least four LBA stations beyond 150 m from the shower axis and where the uncertainties in the amount of circular polarization are less than 0.2. As can be seen from Figs. 1 and 2, there are significant differences between the circular polarization for the thunderstorm events and that for the fair-weather events. First, the circular polarization for the thunderstorm events does not depend on the azimuthal position, φ, of the antenna while for the fair-weather events it is proportional to sin φ. Secondly, the circular polarization for some thunderstorm events changes sign at some distances while the dependence of the circular polarization on distance is almost the same for all fair-weather events as mentioned above. In Fig. 1, it can be seen that there are some thunderstorm events having a very small amount of circular polarization. These events are distinguished from fairweather events by the linear polarization that has been discussed in Ref. [3]. Thirdly, the circular polarization for all fair-weather events is small near the shower axis while it varies from event to event for thunderstorm events. This difference is also shown in Fig. 3, where the amount of circular polarization (jVj=I) within a 30 m radius of the shower axis is given for 884 antennas recording fairweather data and 183 antennas taking thunderstorm data.
We choose the radius of 30 m to concentrate on the near-axis region while also keeping an area large enough to contain a sufficient number of antennas. The uncertainties indicated in Fig. 3 are determined from a Monte Carlo procedure. For 500 trials per antenna the Stokes parameters Q t , U t , and V t are chosen randomly from a Gaussian distribution where the mean and the standard deviation of the distribution correspond to the actual measurement. The Stokes I t of each trial is calculated by using where W is calculated from the actual Stokes parameters measured by the antenna, The spread (standard deviation) of the determined distribution of jV t j=I t gives the uncertainty. Figure 3 shows that the amount of circular polarization near the shower axis is consistently small for fair-weather events, while a large spread is seen for thunderstorm events. In Ref. [8] it was shown that for the fairweather events the measured circular polarization is well understood. The physics of the measured circular polarization of the thunderstorm events is explained in detail in the following section.

III. MODELING
During thunderstorms the emission of radio waves from air showers is affected by atmospheric electric fields [3,14,15]. The atmospheric electric field can be decomposed into two components E ⊥ and E ∥ , which are perpendicular and parallel to the shower axis, respectively. E ∥ increases the number of either electrons or positrons, depending on its orientation, and decreases the other [3,16]. Since the field compensates the energy loss of low-energy electrons, they "live" longer and can thus trail further behind the shower front. As a result, the radiation from these particles does not add coherently in the frequency range 30-80 MHz of the LOFAR LBAs. The transverse component of the field E ⊥ does not change the number of electrons and positrons, but changes the net transverse force acting on the particles [3,16] Hence, the magnitude and the direction of the induced transverse current change according to the net force F ⊥ . Since for the presented data the influence of the transverse component E ⊥ on the radio emission dominates, the parallel component E ∥ is set to 0 in this work. The transverse electric field changes the direction of the transverse current, so it also modifies the polarization of the transverse-current radiation. In thunderclouds, not only the magnitude but also the orientation of electric fields changes with altitude [17]. This causes a change of the transverse current in the thunderclouds and thus the linear polarization changes with time. As explained in the introduction, this results in a changing linear polarization angle over the duration of the pulse giving rise to a large We use a toy model to show the physics of large circular polarization of the pulses in some of the thunderstorm events. We consider the geometry given in Fig. 4 as an example. A vertical air shower passes through two layers where the electric field in each is constant. The fields are such that the net forces are perpendicular to each other and make an angle φ withê v×B as shown in Fig. 4. The induced current in the shower front is proportional to the number of particles in the shower multiplied by the net force acting on them. The induced currents thus have orthogonal directions in the two layers where the peak of the current occurs at height h m , corresponding to X max of the shower, defined as the atmospheric (slant) depth where the number of airshower particles reaches a maximum. For this case we consider the pulses emitted with a central frequency ω when the shower passes through each layer, where η ¼ η a − η b ¼ ωΔt is the phase difference corresponding to an arrival-time difference Δt between the two pulses for an observer. In thunderstorms, the transverse current is generally enhanced by the atmospheric electric field, so its radiation is much larger than the charge-excess emission and thus we ignore the charge-excess contribution. Therefore, neither ε a , ε b nor η depends on the azimuth angle of the antenna position with respect to the shower axis. Since the transverse currents in the two layers are perpendicular to each their pulses are polarized in two perpendicular orientations on the ground. These pulses can be expressed as Substituting these into Eq. (1) we obtain the Stokes parameters For the special case when φ ¼ 0, i.e., the net force in the upper layer is alongê v×B and the one in the lower layer is alongê v×ðv×BÞ (see Fig. 4), the phase shift can be derived from the Stokes V and U parameters We show for this special case how η depends on the distance d from the shower axis for a fixed frequency ω. To simplify the calculation, we assume that at the height of h a the current points inê v×B and emits radiation. After that the shower propagates down with the velocity c, and the current rotates toê v×ðv×BÞ at the height of h b ¼ h a − Δh and radiates another signal. The pulses emitted at different heights move with the reduced velocity v ¼ c=n and thus arrive with a time delay due to the fact that index of refraction n is larger than unity. The signals with a frequency ω which an observer at a distance d from the origin receives (since φ ¼ 0) are ε v×B ¼ A a e iωðt−R a =vÞ ; where are the distances from the observer to the emission points and v is the velocity of the signals. Δh=c accounts for the later arrival of the current at h b . The phase shift of these two signals can be derived from Eq. (7), The phase shiftη is positive, which means that the signal radiated at h a arrives earlier than the one at h b . Forη ¼ 0, the two signals arrive at the observer at the same time. Note that Eq. (9) can only be used in the case where the two emission components are perpendicular to each other and one of the components is alongê v×B . For comparison with the analytic calculation, we simulated three vertical showers with CoREAS [18] that included two-layered electric fields with the boundaries between electric fields at different altitudes h L . The electric field EFIELD option [19] was implemented in CORSIKA [20]. The electric fields in the two layers are such that the net force in the upper layer points inê v×B and the one in the lower layer points inê v×ðv×BÞ , which introduces two perpendicular transverse currents. The upper layer, with strength jE U j¼ 50 kV=m, starts at a height h U ¼ 8 km above the ground and extends down to heights of h L ¼ 4, 3 and 2 km for each simulation. At h L the lower layer starts and the field strength decreases to jE L j¼ 25 kV=m. The shower maximum X max ¼ 580 g=cm 2 is the same in all three simulations, corresponding to h m ≈ 4.6 km, which is in the upper layer.
In order to be compared with the analytic calculation where pulses are assumed to emit a central frequency ω, the phase shift η C from the CoREAS simulations in the narrow frequency band, 60-65 MHz, is derived and displayed in Fig. 5. The phase shiftη derived from Eq. (9) is also shown in Fig. 5 for ω ¼ 65 MHz. To simplify the calculation, the refractive index is kept constant at n ¼ 1.00015. Note that the heights h a and h b in Eq. (9) are the average heights from which the dominant intensity is emitted for the two polarization directions and are thus not equal to h U and h L . In the upper layer, the maximum emission occurs at h a ¼ h m . In the lower layer, the height h b depends on the distance from the observer to the shower axis. At large distances, beyond ∼50 m, the maximum emission arrives from h 0 b ¼ h L − X a =ρ, where the air density ρ is approximately [17] ρðhÞ ¼ 1.208 × 10 −3 expð−h=8.4Þ g=cm 3 and X a is the adapting distance varying with heights (see opening angle corresponding to the distance d, the observer receives the dominant signal from h 0 b ¼ d= ffiffiffiffiffiffiffiffiffiffiffiffi ffi n 2 − 1 p . As seen in Fig. 5, it is the distance at which all three lines coincide. At large distances,η is positive, which means the observer receives the signal radiated at h a first and the one at h b later. At about 50 m, the two signals arrive at the observer at the same time. It can be seen from Fig. 5 that the calculation agrees quite well with the simulations, which demonstrates that the source of the circular polarization is well understood. However, for more general geometries of atmospheric electric fields, the layer heights, field strengths and field orientations can only be found through a numerical optimization procedure.

IV. PROBING THE STRUCTURES OF ATMOSPHERIC ELECTRIC FIELDS
As discussed in the previous section, the circular polarization in thunderstorm events is caused by the variation in the orientation of the atmospheric electric fields. Therefore, using the full set of Stokes parameters, i.e., the combination of intensity, linear polarization, and circular polarization, allows a more accurate determination of the electric fields in the cloud layers where the air shower passes through than when using only intensity information as in Ref. [3]. To provide more insight into this assertion, we discuss in detail the reconstruction of two thunderstorm events which are called in this work event no. 1 and event no. 2.
Fitting thunderstorm events is challenging since the electric fields contain many parameters. Another problem is that since CoREAS is a Monte Carlo simulation, two calculations with similar electric fields can give considerably different results due to shower-to-shower fluctuations even when using the same random seed. Therefore, to determine the electric fields, we first perform a fit using a semianalytic calculation [21] of the radio footprint of air showers based on the current profile. This procedure requires much less CPU time and there are no showerto-shower fluctuations. This allows for a standard steepest descent fitting procedure. Since this method only approximates the structure of the shower front, we use this to get close to the optimal choice after which we use CoREAS for the final calculations. In order to obtain a prediction of the two-dimensional footprints of the four Stokes parameters, we run CoREAS simulations for 160 antennas which form a star-shaped pattern with eight arms as in Ref. [22], and make an interpolation to reconstruct the full profile. The results are filtered in the frequency range of the LOFAR LBAs.
The electric field fields are labeled with indices 1, 2 and 3 where 1 is the top layer. Each layer is defined by the height h above the ground where the electric field starts and the field E ⊥ . Note that our analysis cannot determine the parallel components of the electric fields E ∥ ; therefore we always work in the two-dimensional plane perpendicular tô e v . In this plane, the perpendicular components, E ⊥ , are expressed in two bases. (1) It can be expressed as the field strength jE ⊥ j and the angle α between the net force and e v×B , where the net force is the vectorial sum of the Lorentz force and the electric force given by the electric field. (2) It can also be decomposed into E v×z and E v×ðv×zÞ , the components of E ⊥ alongê v×z andê v×ðv×zÞ , respectively. Hereê z is vertically pointing up. The intensity footprint, the linear polarization footprint and the circular polarization footprint of thunderstorm event no. 1, measured at 12:38:37 UTC, December 30, 2012, are displayed in Fig. 6. The fractions of Stokes parameters are shown in Fig. 7. The intensity footprint (top panel of Fig. 6) of this event shows a bean shape which is also observed in fair-weather events. The differences are that the maximum intensity is not in the v × B-direction as it is in fair-weather events and the linear polarization (middle panel Fig. 6) is not oriented mainly alongê v×B as it is in the fair-weather events. The polarization footprint shows a "wavy" pattern near the shower axis where the polarization is different from the one at the outer antennas. We observe a large fraction of circular polarization in this event, varying as a function of the distance from the antenna to the shower axis. This can be seen in the bottom panel of Fig. 6 and the right panel of Fig. 7. Therefore, for this event, using only the intensity footprint gives incomplete information about the atmospheric electric field. The simplest structure of the electric field which can capture the main features of this event is a three-layer field. The reconstruction is optimal for the values of the parameters given in Table I. The simulation has values of X max ¼ 665 g=cm 2 . The primary energy of the shower is E ¼ 4.7 × 10 16 eV and the zenith angle is θ ¼ 15.5°. Since we do not observe a ringlike intensity pattern, the emission from different layers should not interfere destructively, and thus the fields should not have opposite orientation as taken in Ref. [3]. The change in the orientation of the electric field between the second layer and the third layer, close to the ground, results in a change in the direction of the transverse current and thus gives rise to the rotation of the linear polarization as well as a large amount of circular polarization in the region close to the shower axis. Near the shower axis the radio signal is most sensitive to the later stages of the shower development, while at large distances the currents higher in the atmosphere have more weight. Thus, a much smaller circular polarization component is observed at larger distances. There are some differences between the measured and simulated Stokes parameters seen in Fig. 7 since the three-layered electric field is still an oversimplification of the realistic field. The reduced χ 2 for a joint fit of both the Stokes parameters and the particle data χ 2 =ndf ¼ 4.5, which is large compared to χ 2 =ndf ≈ 1 found in fair-weather showers. However, all the main features are captured.  Figure 8 shows the intensity footprint and the polarization footprint of thunderstorm event no. 2, which was also presented in Ref. [3]. The fractions of Stokes parameters are shown in Fig. 9. This event was measured at 14:28:19 UTC, August 26, 2012. The ringlike structure in the intensity footprint (top panel of Fig. 8) and the overall polarization direction (middle panel of Fig. 8) indicate that at least a twolayered electric field is needed [3], where the electric fields are pointing in opposite directions to introduce a destructive interference between the radiation from the two layers. However, the large amount of circular polarization near the shower axis (see the bottom panel of Fig. 8 and the right panel of Fig. 9) cannot be reproduced by such a field configuration since there is no rotation of the current. The simplest structure of an electric field which can capture the main features of this event is a three-layered field. Table II presents the values of the electric field giving the best reconstruction of this shower. The electric fields obtained here follow the same general structure as presented in our earlier work [3]. Like in Ref. [3] the strength of the fields in the lower layer is about half the strength as in the upper layer with almost opposite orientation. However, in the present, more detailed, analysis an additional layer needed to be introduced which shows that the method used in this work gives more accurate information about the electric fields in thunderstorms. The shower maximum is X max ¼ 628 g=cm 2 . The primary energy of this shower is E ¼ 3.1 × 10 16 eV and the zenith angle is θ ¼ 24.8°. There is also an almost complete reversal of the electric field from the second layer to the third layer which gives rise to the ringlike structure in the intensity footprint and keeps the linear polarization unique. The reduced χ 2 for a joint fit of both the Stokes parameters and the particle data is χ 2 =ndf ¼ 3.5, which is large but reproduces all the main features.
We have checked that the fit quality is sensitive to the heights of the layers on the order of a hundred meters and the orientations of the electric fields at the level of degrees. However, it is not sensitive to heights above 8 km because at that height there are few particles in the shower and thus their contribution to the total radio emission is small. The electric fields shown in Tables I and II only include the components of the true fields perpendicular to the shower axis. The parallel component of the electric fields hardly affects the LBA observations and thus it cannot be determined. In addition, in the frequency domain of the LBAs there is no sensitivity to the component of electric fields in excess of about 50 kV=m, so the strength of the perpendicular component can only be probed up to about this strength (see Ref. [16] for the discussion). To increase the sensitivity, we would need lower-frequency antennas.
However, as explained in the following we have measured large horizontal components of the electric fields along the  shower axis in thunderclouds. A strict vertical electric field can be decomposed into two components, one alongê v and the other one alongê v×ðv×zÞ . Measuring a component in e v×ðv×zÞ (see Tables I and II) could thus be a reflection of a vertical field since the present observations have little sensitivity to anê v component of the electric field. However, a nonzero component in theê v×z direction (see Tables I and II) can never be a projection of a purely vertical electric field, and is thus a genuine signature of a horizontal component. We have confirmed that setting any of the E v×z components to 0 results in poorly reconstructed Stokes parameters. Therefore, it can be concluded that the atmospheric electric field is not fully vertical, but has a significant horizontal component. A three-layered structure and a horizontal component of the electric fields in thunderclouds have also been observed in balloon experiments [2,17,23,24]. The large component of a horizontal electric field at high altitudes can be given by two oppositely charged regions inside a thundercloud. The small horizontal component at low altitudes can be given by the main negativecharge layer of a thundercloud in the center and a local positive-charge region at the bottom of the cloud.

V. CONCLUSION
Air showers measured with the LOFAR LBAs during thunderstorms have generally a much stronger circular polarization component near the shower axis than showers recorded during fair weather. We demonstrate on the bases of a simple model that this is a reflection of the fact that the orientation of atmospheric electric fields changes with height. This gives rise to a rotation in the direction of the transverse current as the air shower proceeds towards the surface of the Earth. This is also confirmed by CoREAS simulations.
Using the full set of the Stokes parameters thus strongly improves the determination of the atmospheric electric fields in thunderclouds. As specific examples we have analyzed two thunderstorm events where we show that the intensity and polarization signature can only be described by a three-layered electric field. Also in balloon measurements, generally three different layers are observed below a height of 8 km. In our analysis, we also determine that the atmospheric electric field has a sizable horizontal component.