key: cord-0434193-p5ie2u5j
authors: Hanson, Jordan C.; Hartig, Raymond
title: Complex Analysis of Askaryan Radiation: A Fully Analytic Model in the Time-Domain
date: 2021-06-01
journal: nan
DOI: nan
sha: f61fa7729235931c8eaa77ccac7b7ddae303d77c
doc_id: 434193
cord_uid: p5ie2u5j

The detection of ultra-high energy (UHE, $>$10 PeV) neutrinos via detectors designed to utilize the Askaryan effect has been a long-time goal of the astroparticle physics community. The Askaryan effect describes radio-frequency (RF) radiation from high-energy cascades. When a UHE neutrino initiates a cascade, cascade properties are imprinted on the radiation. Thus, observed radiation properties must be used to reconstruct the UHE neutrino event. Analytic Askaryan models have three advantages when used for UHE neutrino reconstruction. First, when analytic models are matched to observed data, cascade properties may be calculated directly from single RF waveforms. Second, fully analytic models require no Monte Carlo simulation of cascade particle trajectories, minimizing computational intensity. Third, fully analytic models can be embedded in firmware to enhance the real-time sensitivity of detectors. We derive a fully analytic Askaryan model in the time-domain given the energy and geometry of the UHE neutrino-induced cascade. We compare the fully analytic model to a semi-analytic parameterization used commonly in NuRadioMC, a simulation being used to design the radio component of IceCube-Gen2. We find correlation coefficients greater than 0.95 between the fully analytic and semi-analytic cases, and we find the total power in the signals agree to within 5%.

The flux of neutrinos originating from outside the solar system with energies between [0.01-1] PeV has been measured by the IceCube collaboration [1] . Previous analyses have shown that the discovery of UHE neutrinos (UHEν) will require an upgraded detector design with a larger effective volume because the flux is expected to decrease with energy [2] [3] [4] [5] [6] . Neutrinos with energies above 10 PeV could potentially explain the origin of UHE cosmic rays (UHECR), and provide the chance to study electroweak interactions at record-breaking energies [7, 8] . Utilizing the Askaryan effect greatly expands the effective volume of UHE-ν detector designs, because this effect offers a way to detect UHE-ν with radio pulses that travel more than 1 kilometer in sufficiently RF-transparent media such as Antarctic and Greenlandic ice [9] [10] [11] .

The Askaryan effect occurs within a dense medium of refractive index n when a relativistic particle with v > c/n initiates a high-energy cascade with negative total charge. The net charge in the cascade radiates in the RF bandwidth, and the radiation may be detected by RF channels if the medium is RF-transparent [12] [13] . If the primary particle initiating the cascade is a UHE-ν, the medium must have both sufficient RF transparency and sufficient volume to capture the signals. The Ice-Cube EHE analysis constrains neutrino models predicting a three-flavor neutrino flux of E 2 ν φ ν ≤ 2 × 10 −8 GeV cm −2 s −1 sr −1 between [5 × 10 !5 − 2 × 10 19 ] eV. [4] . Arrays of O(100) in situ detectors that encompass effective areas of ≈ 10 4 m 2 sr per station, spaced by approximately one RF attenuation length, have the potential * Electronic address: jhanson2@whittier.edu to discover a UHE-ν flux below current differential limits. Because ice formations with sufficient volume are abundant in Antarctica and Greenland, a field of prototype Askaryan-class detectors has risen in Antarctica and Greenland. These detectors are designed to discover increasingly constrained UHE-ν flux from astrophysical and cosmogenic sources [5, 6, 14, 15] .

Askaryan radiation was first measured in the laboratory in silica sand, and later ice [16] [17] [18] . The effects that govern the amplitude and phase of the radiation are understood. At RF wavelengths, individual cascade particles radiate coherently, and the overall amplitude scales with the total track length of the excess negative charge. The longitudinal length along the cascade axis determines the RF pulse shape, and the radiation must be observed near the Cherenkov angle. The excess charge profile of a cascade is a function that describes the excesse negative charge versus position along the cascade axis. The radiation procedes from the region near the cascade maximum. At energies far above 10 PeV, however, the charge excess profile can have many local maxima due to the LPM effect [19, 20] . Wavelengths that are shorter than the lateral width of the cascade, perpendicular to the cascade axis, are attenuated. The lateral cascade width has the effect of a low-pass filter on the radiation [20] . Thus, a theoretical foundation has been constructed from a variety of experimental results.

The field of Askaryan-class detectors requires this foundation for at least two reasons. First, the theoretical form of the Askaryan electromagnetic (EM) field is used to optimize RF detector designs. Askaryan models are incorporated into simulations [21] [22] [23] in order to calculate the expected signal in a real detector. Recent Askaryan reconstruction tools for the radio component of IceCube-Gen2 combine machine learning and insights from Askaryan radiation physics [24] [25] [26] . Second, theo-retical models may be used as templates to search large data sets for signal candidates when combined with RF impulse response [5, 27] . At the UHE-scale, the EM fieldstrengths of kiloVolts per meter are expected a few meters from the cascade. However, the detected signal-to-noise ratios (SNRs) at RF channels are expected to be small (SNR ≈ 2) because the typical vertex separation is ≈ 1 km. The amplitude of the EM field decreases with the distance to the vertex (1/r), and there is some alteration of the signal by the ice [9, 28, 29] . The combination of low SNR and RF sampling frequency for the signal bandwidth of [0. GHz implies that RF channels are triggered at high rate by RF thermal noise. The thermal trigger rate causes potential UHE-ν signals to be hidden within millions of thermal trigger waveforms. Thus, template waveform analysis and other digital signal processing approaches must be used to isolate rare signals.

Askaryan models fall into three categories: full MC, semi-analytic, and fully analytic. The original work by E. Zas, F. Halzen, and T. Stanev (ZHS) [13] was a full MC model. The basic properties of the EM radiation from cascades up to 1 PeV were reproduced with threshold energies of 0.5-1 MeV. Sub-threshold particles were assumed to lose energy to processes like the photo-electric effect and not contribute to the EM radiation. Overall signal amplitude was found to scale with energy because energy scales with the total weighted, projected track length. A parameterization for the EM field below 1 GHz was offered, attenuating modes above 1 GHz via a frequency-dependent form factor. The form factor cutoff frequency may be related to the lateral cascade width [20] . The semi-analytic approach was introduced by. J. Alvarez-Muñiz et al (ARVZ) to account for the non-Gaussian fluctuations in the charge excess profile caused by the LPM effect, while keeping the form factor analytic when the EM radiation is observed at the Cherenkov angle [19] . The final EM wave is obtained from the negative derivative of the vector potential once it has been convolved with the charge excess profle, accounting for the viewing angle. The vector potential at the Cherenkov angle is labeled the form factor. Recent work also accounts for differences in fit parameters from simulated EM and hadronic cascades, and other interaction channels, while matching full MC field amplitudes to within a few percent [30] .

Finally, analytic models that reproduce the features of Askaryan radiation from first principles have been introduced. J. Ralston and R. Buniy (RB) introduced a fully analytic model valid for cascades with longitudinal length greater than or less than the length over which EM radiation is coherent, encapsulated by a parameter η [31] . The goal was to create a model valid in the near and far-fields of the EM radiation. Recently, a model was introduced by J. C. Hanson and A. Connolly (JCH+AC) that built upon RB by providing an analytic form factor derived from GEANT4 simulations, and accounted for LPM elongation of the longitudinal cascade width [20] . The form factor and the RB coherence effect formed a set of poles in the complex plane that filter the EM radiation in the same way a multi-pole low-pass filter governs analog signals. The LPM treatment in JCH+AC increased the longitudinal length at energies far above the LPM energy (0.3 PeV in ice), but did not account for the deviation from Gaussian behavior caused by sub-showers in the LPM effect. Both RB and JCH+AC expressed the model generally in the Fourier domain. However, an example EM field was provided in the time-domain if the field is observed at the Cherenkov angle. The theme of this work is to present a fully analytic time-domain model valid for all viewing angles, provided that η < 1.

In Section II, the cascade geometry, units, and relevant vocabulary are defined. In Section III, we present the derivation of the form factor from first principles. Appendix A of the JCH+AC publication contains the detailed comparison between GEANT4 cascade simulations and theoretical predictions for the lateral charge distribution. In Section IV, a new derivation of the Askaryan EM field observed at the Cherenkov angle (on-cone) is presented that accounts for the relative strength of one complex pole from coherence effects and two from the form factor. In Section V, a new derivation of the Askaryan EM field not observed at the Cherenkov angle (off-cone) is presented that accounts for the longitudinal cascade width, viewing angle, and form factor. In Section VI, analytical results are compared to published Askaryan EM fields generated by the semi-analytic parameterizations in NuRadioMC [23] at 10 PeV and 100 PeV. At this energy, the LPM effect only has a small influence on the cascade development. In Section VII, the results are summarized and several applications of the new analytic model are described.

The coordinate system of the Askaryan radiation from a vector current density J is shown in Fig. 1 (a)-(b). Primed coordinates refer to the cylindrical coordinate system of the high-energy cascade. The zenith or viewing angle is measured with respect to the longitudinal axis (z ). The separation between the origin and observer is r = | x − x |, in ther direction in spherical observer coordinates. The origin is located along the longitudinal axis where the cascade has the highest instantaneous charge density (ICD). The cascade is assumed to generate an ICD with cylindrical symmetry and therefore no dependence on φ . This assumption is based on the large number of cascade particles, and the conservation of momentum. The lateral extent of the ICD is along the lateral axis (ρ ). The viewing angle in Fig. 1 (a)-(b) is θ in spherical coordinates, and the Cherenkov angle occurs when θ satisfies cos(θ C ) = 1/n ice . The refractive index for solid ice is n ice = 1.78 ± 0.003 [32] .

In Fig. 1 (c) , an archetypal cascade profile is shown in the form of an excess charge function n(z ) with char- 1: (a) Side view of the coordinate systems used in the analysis. Coordinates without primes are from a spherical system and refer to the observer. Primed coordinates refer to the reference frame of the vector current density, J(ρ , z ), and are from a cylindrical system. (b) Front view of the coordinate system. The instantaneous charge density (ICD) is assumed to have cylindrical symmetry with no φ -dependence. (c) The function n(z ) describes the total cascade excess charge, and it has a characteristic width a. At any given moment, the ICD has a width much smaller than a [20] .

acteristic longitudinal length a. The individual ICDs represent the excess charge density in moments in time, whereas n(z ) refers to the excess charge as a function of z = ct . Approximating the central portion of n(z ) as a Gaussian distribution N (µ, σ) corresponds to setting a = 2σ. Askaryan radiation occurs because of the excess charge from cascade electrons (≈ 20%), an effect verified by ZHS and later with GEANT4 [13, 20, 33] . Cascades may be characterized as electromagnetic, initiated by charged outgoing leptons from the UHE-ν interaction, or hadronic, initiated by the nucleus with which the UHE-ν interacts. In fully analytic models, cascade profiles are given by the Greisen distribution (electromagnetic case) and the Gaisser-Hillas distribution (hadronic case). This is also true in MC codes like AraSim that utilize the semi-analytic ARVZ parameterization [11] . The units of the electromagnetic field in the Fourier domain are V/m/Hz, often converted in the literature to V/m/MHz. To make the distance-dependence explicit, both sides of field equations are multiplied by r, as in r E = ..., making the units V/Hz. Throughout this work, an overal field normalization constant E 0 is used in the time domain, and it is directly proportional to the cascade energy E C , as in most Askaryan models that followed from the original full MC results [13] . Thus, E 0 may be linearly scaled with energy, provided the longitudinal length a and the maximum number of cascade particles are derived consistently from the appropriate cascade distribution. When the field amplitude is proportional to E 0 times a characteristic frequency-squared (on-cone), the units of E 0 must be V/Hz 2 . When the field amplitude is proportional to E 0 times a characteristic frequency divided by a characteristic pulse width (off-cone), the units of E 0 remain V/Hz 2 .

It will be shown in Section III B that the longitudi-nal length a in both the electromagnetic and hadronic cases demonstrates a weak dependence on the ratio of the initial cascade energy to the cascade critical energy: a ∼ ln (E C /E crit ). Specifically for the Greisen distribution, it can be shown that if n max = n(z max ), where z max = ln (E C /E crit ), then n max a ∼ E C /E crit . Thus, the area under the curve n(z ) scales as E C . Among the RB results is the idea that the Askaryan radiation amplitude is proportional to n max a. Although the amplitude develops over a cascade of length a, the RB field equations remain valid for a coherence length ∆z coh in which r(t) is stable to first order relative to a wavelength. The η parameter is the square of the ratio of a to ∆z coh :

The dominant Askaryan radiation occurs when |z | ∆z . In far-field, η < 1. The calculations, however, are valid for any η, rather than only the in the far-field (η 1, and kr 1, r a). If a ∆z coh , then the fields have spherical symmetry, and the limit kr 1 corresponds to the far-field. Conversely, if a ∆z coh , and kr 1, then η 1 and the radiation is observed in the near-field. In the first JCH+AC model, a limiting frequency ω C (Equation 2) was shown to filter the Askaryan radiation. The effect of ω C on the Askaryan E-field is described in Section IV.

The Askaryan radiation is primarily polarized along the unit vectorθ, although there is a small amount of amplitude alongr [20] . The wavevector associated with propagation along r and polarized linearly alonĝ θ is k = (2π)/(nλ), where n is the index of refraction. A 3D wavevector is defined by RB, equivalent to q = nk(1, ρ/R). Let the vector current density be equal to the charge density times the velocity of the ICD: J(t, x ) = ρ(z − vt, ρ ) v. Further, the charge density is factored into charge times ICD as a number density: ρ(z − vt, ρ ) = n(z )f (z − vt, ρ ). The form factor that filters theθ-polarized radiation is the three-dimensional spatial Fourier transform of f (z − vt, ρ ):

The result for F was calculated by JCH+AC [20] , and more detail on that derivation is provided in Section III A, along with the energy-dependence on the longitudinal length a in Section III B. One main result shown by JCH+AC is that the effect of F on the radiation is that of a two-pole low-pass filter if σ < 1, where σ is the ratio of angular frequency to the cutoff frequency introduced by the form factor:

Armed with the form factor, the longitudinal length parameter, and the RB field equations, the main electromagnetic field may be assmbled according to the following form:

The factor E 0 represents the overall amplitude scale, and it is proportional to cascade energy. The factor ω is the angular frequency. The variable ψ is ψ = −i exp(ikr) sin θ. The function E(ω, θ) contains the vector and complex pole structure of the field (see [20] ). As presented above, F (ω, θ) is the form factor. Thus, the model represented by Equation 5 is an all-θ, all-ω model. Equation 5 is valid at all frequencies not attenuated by F (ω, θ), and all viewing angles not attenuated by E(ω, θ). The goal of the following sections is to build an all-θ, all-t model in the time-domain, derived from Equation 5.

In order to arrive at the main electromagnetic field in the time-domain, the individual pieces of Equation 5 must first be assembled. The first piece will be the form factor F that accounts for the 3D ICD, followed by some remarks about the energy-dependence of the longitudinal length parameter a.

Let an observer be a distance r from a charge distribution f (z , ρ ) in a primed cylindrical coordinate system x = (z , ρ ). Let the scalar wavevector in a material with index of refraction n be k = 2π/(nλ), and the angular frequency be ω = ck. Recall the definition of q from Section II: q = nk(1, ρ/R). As presented in Section II the form factor for the cascade ICD is the 3D Fourier transform of f (z , ρ ).

The goal is to evaluate F in the Fourier domain for an ICD definition informed by cascade simulations. Simulations of the cascade induced by the neutrino indicate a thin wave of charge in z spread uniformly in φ , that decreases exponentially in ρ [20] . Thus, The thin wave of charge moves at speed c, and the ρ 0 parameter has units of inverse length. The ICD function f (z , ρ ) is a number density, and integrating over ρ dρ dφ , with z = ct , shows that f is normalized to 1. The full derivation of F given f (z , ρ ) is given in Appendix A. The final result is a simple analytic formula:

According to Equation 8 , the form factor is a low-pass filter with the cutoff-frequency ω CF introduced in Section II: σ = ω/ω CF :

F ≈ 1

The definition ω 0 = 2/3 ω CF has been used in the final step. The two-pole approximation resembles the original ZHS parameterization that matched full MC results below 1 GHz (see Equation 20 of [13] ).

In Section VI below, the parameter l = 1/( √ 2πρ 0 ) will be derived from best-fit ω 0 values. The parameter l is the decay constant of the exponential in the ICD (Equation 7). The l-parameter can vary because the best-fit ω 0 varies, and the two numbers are connected by Equation 4 . It is important to note that l = R M , where R M is the constant Molière radius.

The Molière radius is defined as the lateral radius which forms a cylinder containing 90% of the energy deposition of a cascade in dense media. Let X 0 be the electromagnetic radiation length in ice, and let Z be the atomic number. A common formula for the Molière radius, R M , is

Another useful formula relates X 0 to Z and the mass number A:

Using Z = 8 and A = 16 for oxygen: X 0 = 34.46 g cm −2 . Converting R M to units of cm, with a density of 0.917 g cm −3 , gives R M = 9.16 cm. The measured value for water is 9.77 cm. Although it is tempting to compare l to R M , these parameters have different definitions. Knowing that l is related to ω 0 , l may be estimated as λ/2 at the cutoff-frequency in ice. At 3 GHz in ice, λ/2 ≈ 2.8 cm, and at 1 GHz in ice, λ/2 ≈ 8.4 cm. Although the results are at the same order of magnitude as R M , it will be shown that there are three effects limiting the high-frequency spectrum of the radiation: ω 0 , ω C , and the viewing angle. Thus, l can be less than R M and the radiation can still be limited to 1 GHz.

The next piece required in the assembly of the main electromagnetic field is the energy-dependence of the overall amplitude, and the energy dependence of the longitudinal length parameter, a. What follows are two separate discussions, one for electromagnetic cascades, and one for hadronic cascades.

The number of charged particles versus distance n(z ) in an electromagnetic cascade taking place in a dense medium with initial cascade energy E C , critical energy E crit , normalization parameter n 0 , and age s is [20] 

The age is given by

The function n(z ) is known as the Greisen distribution. In this case, the longitudinal distance variable z is in radiation lengths, and n(z ) is maximized at

The calculation of the energy dependence of a is broken into four steps: normalization of n(z ) to be a fraction of the maximum, n max , conversion of n(z ) to n(s), determination of the width of n(s), and conversion of the width in s to a width in z .

Substituting s = 1 and z = ln(x) into n(z ) yields n max :

Dividing n(z ) by n max yields

The function f (z ) is the normalized Greisen distribution. Solving the age equation for z :

From this point forward, since f (z ) is being studied near s = 1, we use ln(s) ≈ s − 1. Substituting back into Equation 17 :

To calculate the Gaussian width of f (s), consider the following procedure. Suppose there is a function f (y) such that ln(f (y)) is equal to a negatively oriented quadratic: ln(f (y)) = −(1/2)y 2 /σ 2 . Solving for y yields y ± = ± √ −2 ln f σ. Assuming f (y) behaves like a Gaussian near the global maximum, the width is ∆y = y + − y − = 2 √ −2 ln f σ and ln f = −0.5 precisely when it should: y = σ. In a similar fashion, let q = ln(xf )/ ln(x) in Eq. 19. This leaves

Letting b = (q + 5)/3,

The longitudinal length ∆z is required. The conversion of Equation 23 can be made using Equation 18 . The result is

The expression for ∆z may be simplified using the definitions of s ± , q, and b:

Noting that q = 1+ln(f )/ ln(x), and that the quadratic factors:

The term involving (ln(f )) 2 is small, positive, and negligible compared to the second term in that square root. Dropping it and setting ∆z = a yields the main result:

Recall that f < 1, so ln(f ) < 0 and a is real-valued. Equation 29 is in radiation lengths. In solid ice the density is ρ ice = 0.917 g cm −3 , and the electromagnetic radiation length is z 0 = 36.08 g cm −2 [20] . Converting to distance gives

Note that a ∝ ln(x). Since n max ∝ x/ ln(x), the product n max a, is an approximation of the area under n(z ), and it goes as the energy x = E C /E crit . This was noted by RB [31] , who placed n max a in the electromagnetic field normalization rather than E C . Allowing f to drop to 40% of its peak value, and using E crit ≈ 10 8 eV, gives a ≈ 4 m for E C = 10 16 eV. It will be shown in Section VI that 4 meters is a good number for a 10 16 eV cascade derived from full MC simulations.

The Gaisser-Hillas distribution describes hadronic cosmic-ray air showers, but has also been applied to hadronic cascades in dense media [22] . The original function reads

The variables are defined as follows: N max is the instantaneous maximum number of particles in the cascade, z is the longitudinal distance in radiation lengths, z 0 is the initial starting point, λ is the interaction length, and z max is the location of N max .

The calculation of the energy dependence of a is broken into three steps: the normalization of n(z ) to be a fraction of n max with convenient parameters, approximation the distribution near z max , and the longitudinal width determination that contains a pre-defined fraction of the normalized n max . Let x = (z − z 0 )/λ), and w = (z max − z 0 )/λ. One useful (normalized) parameterization is

Taking the logarithm of both sides of Eq. 32 and then dividing both sides by w gives:

Near the shower maximum, x/w ≈ 1. Let = x/w, and expand the right-hand side in a series about = 1 to second order, noting that the first order vanishes:

Notice that the logarithm of the normalized distribution is again a negative quadratic, indicating Gaussian behavior near z max . Let b = ln(f )/w, and solve the quadratic equation for :

As with the electromagnetic case, let ∆ = + − − , where + and − are the two solutions to the quadratic. The result is

Assuming constant density, the starting point of the cascade z 0 is arbitrary, so z 0 may be set to zero. Also, ∆x = w∆ , so

If z 0 = 0, w = z max /λ. Finally, knowing that ∆z = λ∆x, the result becomes

As with the electromagnetic case, the longitudinal length parameter a goes as √ z max .

Theθ-component of the electromagnetic field at θ = θ C will now be built in the time-domain, using the general RB field equations in the Fourier domain, along with the form factor and longitudinal length parameter. Beginning with Eq. 45 JCH+AC [20] , assuming Equation 8 for the form factor, with σ = ω/ω CF and η = ω/ω CF , and letting E 0 be proportional to cascade energy E C :

More detail is provided in Appendix B. Let the retarded time be t r = t − r/c, and let ω 0 = 2 3 ω CF and E 0 = E 0 sin θ C . The inverse Fourier transform of Eq. 39 is (see Appendix B)

In Equation 40, the derivative with respect to the retarded time d/dt r is introduced to remove a factor of (−iω) from the numerator. Further, the three simple poles arise from approximation in the denominator. The poles are located on the imaginary axis, above and below the origin. Thus, region of convergence of Eq. 40 depends on the sign of the retarded time. The details of each case are provided in Appendix B. The results, including the effect of each pole and each sign of t r , become

Equation 41 represents the time-domain solution for the on-coneθ-component of the Askaryan electric field. Table I Table I are shown in Section VI.

Notice that the overall field amplitude scales with energy (E 0 ∝ E C ∼ n max a) and with the second power of cutoff-frequency. The units of E 0 are therefore V Hz 2 . The amplitude is asymmetric, and the the parameter influences the asymmetry. The parameter was studied in JCH+AC in detail. For example, Fig. 10 of [20] shows that ≈ [0.1 − 1] for √ 2πρ 0 ≈ 20 m −1 and a ≈ 4 m. These values are typical of 10 PeV cascades (see Section. VI). The expression for is the product of the ratio of the lateral to longitudinal length, and the longitudinal length to the vertex distance, making it a physical parameter connecting the event geometry to the cacscade shape. Figure 3 displays normalized examples of Equation 41 for different values of ω 0 , ω C , and . The values shown are motivated by comparisons to semi-analytic parameterizations in Sec. VI.

Among the results of JCH+AC is the Gaussian width of the radiation in the Fourier domain: σ ν , where ν represents the frequency. Generally speaking, the product of the width in the Fourier domain with the width in the time domain must be greater than or equal to some constant. Gaussian functions minimize this uncertainty principle. The following procedure is used to compute the width σ t of the on-cone field. First, the branches corresponding to t r < 0 and t r > 0 are each treated as probability distribution functions, and normalized. Next the average positive and negative retarded times,t r,+ and t r,− , are computed: Subtracting the two averages yields σ t , which simplifies to

The result has the correct units and the limiting cases are sensible. Suppose → 0 (ω C ω 0 ), then σ t → 2/ω 0 , which is expected from observing Equation 41 if the ω C exponential disappears. If = 1 (ω C = ω 0 ), then σ t = 3/ω 0 . That is, the pulse is wider if there is more than one relevant cutoff frequency. The result also reveals that as the pulse asymmetry grows with , so does the pulse width.

Equation 36 of JCH+AC [20] provides σ ν :

Expanding to first order in ∆ cos(θ) = cos(θ)−cos(θ C ),

From previous definitions, let ω −1 C = a 2 sin 2 (θ C )/(rc), and ω −1

(47) The limits r a and ka 1 cause η → 0. Assume these limits to be true, along with ∆θ → 0. The result is

Therefore, in order to satisfy σ ν σ t > 1/(2π), Although a/r 1 and l/a 1, as long as these expressions do not approach zero as fast as ∆θ → 0 in Equation 49, the uncertainty principle holds. Yet these are exactly the conditions of the problem. The observer distance r is assumed to be in the far-field and much larger than the longitudinal length, and at cascade energies at or above 10 PeV, the longitudinal length is computed to be much larger than the lateral ICD width. Thus, the result is that σ ν σ t > 1/(2π).

Turning to the case for which θ = θ C , theθ-component of the electromagnetic field will now be built in the timedomain. The RB field equations for theθ andr components are summarized in both RB and JCH+AC [20, 31] . Recall the general form of the electromagnetic field, given in Equation 5: 

The first task is to simplify E(ω, θ) before taking the inverse Fourier transform. The full definition of E is given in Appendix C. The simplification inolves expanding E(ω, θ) in a Taylor series such that u = 1 − iη ≈ 1, restricting η < 1. The restriction η < 1 may be satisfied by either a sin θ/r 1, or ka sin θ 1, or both. Once E(ω, θ) is simplified, the inverse Fourier transform of Equation 50 may be evaluated to produce the result. Table II contains useful variable definitions, Table III contains useful function definitions, and Table IV contains special cases of the functions in Table III .

The original form of E(η, θ) is shown in Appendix B. Changing variables to u and x (Tab. II) and using the 

The details of the expansion are shown in Appendix C. The result is

The inverse Fourier transform of theθ-component gives the time-domain results, after including the expanded E(u, x):

Intriguingly, the result is proportional to the linebroadening function, H (DLMF 7.19, [34] ) common to spectroscopy applications. There are three terms in Eq. 53. Two terms ultimately vanish, being integrals over odd integrands (see Appendix C). The integral that remains contains H, with ω 1 = t r /(2p):

Though the line-broadening function cannot be expressed analytically, there are examples of polynomial expansions [35] . Note that, for situations relevant to the current problem, ω > ω 1 . Requiring that ω > ω 1 amounts to a restriction between ∆θ and |t r |:

For example, taking a = 5.0 meters, a cutoff frequency of ν = 1 GHz, and letting ∆θ be as small as 1.5 • allows |t r | 0.9 ns. For ∆θ = 3.0 • , |t r | 3.6 ns. The field quickly approaches zero outside this window. Hereafter, this step will be called the symmetric approximation. The restriction on ∆θ is formalized in Sec. V B. Using the results of Sec. V A, one can show that |t r |/σ t < ωσ t , or that |t r |σ ω < ω/σ ω . Solving I 0 clears the way for the final result (see Appendix C): Figure 5 contains contours of the same results as in Fig. 4 . Table I Table I are shown in Section VI.

As in the on-cone result, the overall field amplitude scales with energy (n max a). However, the amplitude scales with the first power of cutoff-frequency. The units of E 0 remain V Hz 2 because of the factors of t r and p, which have units of time and time-squared, respectively. The argument of the complementary error function, √ pω 0 , is unitless. This factor is strictly positive, so the range of the complementary error function is (0, 1). The factor √ pω 0 cannot be zero without setting θ = θ C , or setting ω CF = 0. Both cases are not allowed. Equation  57 represents the off-cone (θ = θ C ) solution, so p = 0. Another possibility is that p = 0 if a = 0, but this implies E 0 = 0. Further, ω CF = 0 would require cascade particles to have no lateral momentum, which is not physical. Therefore, 0 < erfc( √ pω 0 ) < 1.

As in Section IV A, the uncertainty principle should be checked. Equation 57 is an anti-symmetric Gaussian function with pulse width σ t = √ 2p. Let ∆ cos θ = (cos θ − cos θ C ). Using Table II , the expression √ 2p evaluates to

Recall that σ ν is given by

The uncertainty product is

Expanding near u ≈ 1 implies that η → 0. This condition corresponds to observations in the far-field. Thus, the model obeys the uncertainty principle. For the special case η = 0, the model minimizes the uncertainty product.

The form of Eq. 57, and the restriction between ∆θ and |t r | from the symmetric approximation suggests the limit ∆θ = θ − θ C must be examined carefully. Since p ∝ (cos θ − cos θ C ) 2 , probing the model near θ = θ C is equivalent to taking the limit that p → 0. Intriguingly, the p −1 -dependence in the field does not lead to a divergence. As the field grows in amplitude from p −1 as p → 0, the field width, √ 2p, approaches zero. Thus, the power in the field does not diverge.

Equations 44 and 58 contain the pulse widths of the on-cone and off-cone fields, respectively. Power in the offcone case is limited by the pulse width √ 2p, and the observed power should decrease in coherence as ∆θ and the pulse width both increase. Thus, a reasonable constraint on ∆θ min is given by setting the minimum off-cone pulse width equal to the on-cone pulse width. Setting Eqs. 44 and 58 equal:

Expanding the expression for p near θ = θ C , and evaluating the square root leads to

Using = ω 0 /ω C , and letting k 0 = ω 0 /c, the formula may be rearranged:

Squaring both sides, and then dividing both sides by r yields

The quantity in parentheses is η, with ω = ω 0 . Since η = ω/ω C , setting ω = ω 0 means η = . Solving for ∆θ min gives

Assuming ≈ 1, f 0 ≈ 1 GHz, n = 1.78 for solid ice, and c = 0.3 m ns −1 (see Sec. VI A), k 0 ≈ 35 m −1 . Taking r = 1000 m, ∆θ min ≈ 1 • . Simple rules to remember for the application of the off-cone field are:

The fully analytic model will now be compared to the ARZ semi-analytic parameterization used in NuRa-dioMC [30] to predict signals in IceCube-Gen2 Radio. To provide concrete comparisons, a small set of waveforms were generated with NuRadioMC, for both electromagnetic (EM) and hadronic (HAD) cascades. The EM cascades have an energy of E C = 10 16 eV, while the HAD cascades have an energy of E C = 10 17 eV. The reason the EM energy is lower is to avoid the LPM effect, which changes the shape of n(z ), and requires more careful modeling [20] .

The comparison involves three stages. First, waveforms and a-values are generated for each cascade type, energy, and angle: θ = θ C + 3.0 • , and θ = θ C . Second, Equations 41 and 57 are tuned to match the waveforms. In each fit, the Pearson correlation coefficient (ρ) between MC waveforms and Equations 41 or 57 is maximized, and the sum-squared of amplitude differences ((∆E) 2 ) is minimized. Finally, best-fit parameters are tabulated and compared with MC values from NuRadioMC.

Two remarks are important regarding the fit criteria. First, the Pearson correlation coefficient is not sensitive to differences in amplitude. The coefficient is normalized by the standard deviation of each signal, and the signals have no DC components:

Thus, the parameters that control correlation are those that scale t r . Second, the parameters that control (∆E) 2 are those that scale the waveform amplitude. Note that

The (∆E) 2 minimization resembles χ 2 minimization, but lacks weighting in the denominator from statistical error in the fully analytic model. No statistical error has been associated with Equations 41 and 57. The utility of minimizing (∆E) 2 is that the final difference in power between model and data is minimized. samples. Let f C = ω C /(2π) and f 0 = ω 0 /(2π). The frequencies f C and f 0 were varied from [0.6 -6.0] GHz. The parameter E 0 was varied from [0.05 -5.0] V GHz −2 . In a simple 2-level for-loop, the Pearson correlation coefficient ρ was maximized by varying f 0 and f C . Next, the sum of the squared amplitude differences (∆E) 2 was minimized by varying E 0 , while holding f 0 and f C fixed. Several other schemes were studied, including a 3-level for-loop, but the two-stage process produced the best results. The results are shown in Fig. 6 .

Typically, maximizing ρ corresponds to minimizing (∆E) 2 , but the waveforms fluctuate in shape and amplitude. In Fig. 7 , (∆E) 2 is graphed versus ρ for one event. Typical best-fit ρ-values are ≈ 0.97 for this set, corresponding to best-fit (∆E) 2 values of ≈ 7%. Contours of ρ > 0.95 for f 0 versus f C are shown in Fig. 6 (Left column). The crosses represent the best-fit location. The dashed gray line at y = x corresponds to f 0 /f C = = 1. Although Equation 41 contains an expansion to first order in , the motivation for the expansion is to cause the formula to resemble the derivative of the vector potential presented by ARZ [30] . The only real restriction is that = 2, due to the singularity. Thus, the best-fit results avoid the solid black lines ( = 2) in Fig. 6 , but find -values large enough to account for pulse symmetry. The best-fit waveforms are shown in Fig. 6 (Right column), where the gray curves correspond to ARZ from NuRadioMC and the black curves represent Equation 41. Table V contains the best-fit results for the Eq. 41 parameters, along with best-fit ρ-values and (∆E) 2 -values. The typical fractional power difference is ≈ 7 %, and correlation coefficients are ≈ ρ = 0.97. For every waveform, f 0 > f C , implying > 1. This result will shift to < 1 below for the hadronic events. To estimate the acceptable range for f 0 and f C , the horizontal and vertical distances from the crosses to the ρ > 0.95 contour is assigned to the errors in Tab. V. Though this error definition leads to large fractional errors in the frequencies and a, the aerrors encompass the a-values from NuRadioMC. Though outside the scope of this work, future work would include a full-scale NuRadioMC run with thousands of waveform comparisons. The region in [f 0 , f C ] space for which UHEν signals are expected for IceCube-Gen2 Radio would be established, and man-made backgrounds that mimic UHE-ν could be rejected based on the best-fit location in the space. Finally, the connection between and cascade type (hadronic or electromagnetic) would be probed.

Hadronic case. Using the same procedure as the electromagnetic case, NuRadioMC was used to generate six hadronic cascades at 100 PeV for comparison to Eq. 41. The energy was increased to show that the waveforms and model scale with energy. The hadronic cascades are not morphed my the LPM effect, so Eq. 41 remains valid. The same fit procedure as the electromagnetic case was used, with two exceptions. First, the range for E 0 was increased by a factor of 10, because the cascade energy is E C = 100 PeV. Second, the correlation contours represent ρ = 0.985. In a simple 2-level for-loop, the Pear- son correlation coefficient ρ was maximized, fixing the complex poles. Next, the sum of the squared amplitude differences (∆E) 2 was minimized with the best-fit pole frequencies.

The main results are shown in Figure 8 . The results shown in Figure 8 demonstrate that modeling hadronic cascades at θ = θ C is similar to the electromagnetic case, except that the contours indicate best-fits with < 1. The amplitudes of the best-fit analytic functions in Figure 6 are smaller by a factor of 10 due to the smaller E C . Recall that the amplitude of Eq. 41 goes as E 0 ω 2 0 , so those parameters will combine to produce the factor of 10 difference. The gray contours in Fig. 8 cor- respond to ρ = 0.985, because the fits are slightly better compared to the electromagentic case. Table VI contains the best-fit parameters corresponding to Figure 8 , along with results for ρ and (∆E) 2 . The typical power difference (∆E) 2 has decreased significantly with respect to the electromagnetic case, despite using identical fit procedures. The ρ-values all exceed 0.985, and the (∆E) 2 results are typically below 2 percent. Intriguingly, < 1 means higher f C values, which in turn yields systematically low a-values, despite the increased energy. Reconstructed a-values are still within a factor of 2 of those derived from NuRadioMC. The parameter controls pulse asymmetry, and the hadronic pulses with < 1 appear to be more asymmetric than the electromagnetic. Despite systematic offsets, the best-fit a and the NuRadioMC a value are tightly correlated (see Fig. 11 below) .

Electromagnetic case. Six different electromagnetic cascades and the corresponding Askaryan fields were generated with NuRadioMC for comparison to Equation 57. The events had θ = θ C + 3.0 • , but r and E C remained unchanged relative to Sec. VI A. The fit procedure is similar to Section VI A. One difference is that ω 0 only changes the waveform amplitude, along with E 0 . The pulse width σ t = √ 2p connects the longitudinal length a and the viewing angle with respect to the Cherenkov angle |∆ cos θ| ≈ sin θ C (θ − θ C ).

The fit procedure was done in two stages. First, θvalues and a-values were scanned from [θ C + 1.5 • , θ C + 10.0 • ] and [0.1, 10] meters, respectively, to determine the best-fit ρ. Once the best-fit values for a and θ are determined, (∆E) 2 is minimized by varying f 0 = ω 0 /(2π) and E 0 from [0.3, 3.0] GHz and [0.1, 2.0] V GHz −2 , respectively. First, θ-values and a-values were scanned in a 2-level for loop. Once the best-fit values for a and θ were determined, (∆E) 2 was minimized by varying f 0 and E 0 in a second 2-level for loop. The main results are shown in Figure 9 .

In Figure 9 (left column), the best-fit a-values and θvalues are marked with a cross. The circles represent the MC values from NuRadioMC. Circles and crosses lie on the dashed lines, because an uncertainty principle connects a-values to θ-values (see Section V A). Specifically, Equation 58 may be used to show, to first-order in ∆θ = θ − θ C :

Treating p as a constant derived from the waveform related to pulse width (σ t = √ 2p) means that the product of a and ∆θ is a constant. The parameters are therefore inversely proportional: a = const/∆θ. Further, the ρ > 0.95 contour is centered along the dashed lines, indicating all good fits respect the uncertainty principle. Thus, measurement of the true pulse width (beyond that introduced by RF detection components) represents a strong constraint on the event geometry. The best-fit results happen to fall above and below the true values along the dashed line, indicating no systematic bias. The best-fit waveforms are shown in Figure 9 (right column). Typical correlation coefficients exceed ρ = 0.98 and typical fractional power differences are less than 3%. Table  VII contains the results from the fit procedure.

The results in Table VII demonstrate that Equation  57 is in excellent agreement with the semi-analytic parameterization. On average, the fits encompass the MC values of θ and a, and the results for f 0 and E 0 values match expectations. The results for l, the lateral width parameter, are derived from the f 0 results. Cascades with l-values of a few cm should have cutoff-frequencies in the ≈ 1 GHz range. The ρ-values are comparable to the hadronic on-cone data (ρ ≈ 0.98), indicating that the symmetric approximation made in Section V does not incur much penalty. The fractional power differences (∆E) 2 are ≈ 2 %, most likely due to the amplitude symmetry of the fully analytic model.

Hadronic case. Six different hadronic events and the corresponding Askaryan fields were generated for comparison to Equation 57. The fit procedure is the same as the electromagnetic case, except that the range for E 0 is expanded to [1.0, 20.0] V GHz −2 , because the cascade energy is E C = 100 PeV. The main results are shown in Figure 10 .

As with the electromagnetic case, ρ is maximized and (∆E) 2 is minimized. Table VIII contains From left to right, the viewing angle, longitudinal length parameter, form-factor cutoff frequency, the energy-scaling normalization, the lateral width of the cascade l = 1/( √ 2πρ0), the best-fit correlation coefficient, and the relative power difference between NuRadioMC semi-analytic parameterization and the fully analytic model. The parameter means and errors in the mean are quoted in the bottom two rows. ρ ≈ 0.98 and (∆E) 2 ≈ 5 % were found. Similar to the results shown in Table VII, the results in Table VIII are in agreement with the MC values from NuRadioMC. Note that the f 0 and l values are statistically equivalent between Tables VII and VIII. The similarity leaves the amplitude scaling to the E 0 parameter. The E 0 -values match expectations for 100 PeV cascacdes, because they are a factor of 10 higher than those of the 10 PeV electromagnetic case.

One interesting question is whether the analytic model has potential to reconstruct cascade parameters from the waveform shape. What follows is a simple reconstruction exercise. For the first exercise, θ = θ C + 3.0 • is assumed fixed by another part of the analysis. For example, θ could be determined by measuring the cutoff-frequency in the Fourier domain below 1 GHz [20] . Fixing θ means a is the only parameter varied. Scanning Equation 57 over all NuRadioMC waveforms at fixed θ = θ C + 3.0 • yields Figure 11 , in which the best-fit a-value is graphed versus the MC a-value. The a-errors from the digitized MC cascades profiles, and the a-errors from waveforms, are taken to be ±10 cm (± two ∆a step-sizes). A nonlinear least-squares (NLLS) algorithm was applied to fit a function to the data, and a linear correlation appeared. The linear function is a good fit, and the correlation be- form-factor cutoff frequency, the energy-scaling normalization, the lateral width of the cascade l = 1/( √ 2πρ0), the bestfit correlation coefficient, and the relative power difference between NuRadioMC semi-analytic parameterization and the fully analytic model. The parameter means and errors in the mean are quoted in the bottom two rows.

tween waveform a-values and MC a-values is 0.97.

The results in Figure 11 imply a reconstruction technique for ln(E C ) using the formulas found in Section III B. Consider the relationship between a and ln(x), where x = E C /E crit : a = c 1 ln(x). The fractional error in ln(x) is proportional to the fractional error in a:

To advance the field of UHE-ν detection, a fully analytic Askaryan model in the time-domain has been presented, that accounts for the UHE-ν cascade energy, and lateral and longitudinal geometry. The fully analytic model was compared to published results generated with a semi-analytic parameterization used commonly in Nu-RadioMC, a simulation being used to design IceCube-Gen2. Pearson correlation coefficients between the fully analytic and semi-analytic paremeterizations were found to be greater than 0.95, and typical fractional differences in total power were found to be ≈ 5%. What follows is a brief summary of new results, and a list of potential uses for a fully analytic model in UHE-ν detection. 

The main results are summarized in Table IX . Beginning with the analytic form factor, it was shown that the lateral component of the cascade ICD leads to the effect of a two-pole, low-pass filter on the radiated Askaryan field. The analytic form factor was first presented in [20] , however this work has introduced the first complete study of the effect of its complex poles on the Askaryan radiation. The longitudinal factor a was introduced by Ralston and Buniy [31] , who described the connection to cascade energy. In this work, a technique for cascade energy reconstruction was presented, based on constraining a-values by fitting analytic functions to observed waveforms. The connection was made for both electromagnetic cascades and hadronic cascades via common parameterizations of each cascade type.

The on-cone field equations were derived, similar to the final section of [20] . This work represents the first time the two pole frequencies f 0 and f C and the asymmetry parameter have been used to characterize the time-domain field equations such that they match semianalytic parameterizations. The parameter reveals a potential cascade type identification scheme. Finally, the on-cone field equations were shown to satisfy the uncertainty principle for Fourier transform pairs. The off-cone field equations were derived as well, and were shown to be in excellent agreement with semi-analytic parameterizations. Cascade parameters like a and the viewing angle θ were fit directly from the waveforms. Finally, matching off-cone field equations to data shows promise as an energy reconstruction.

To obtain the off-cone field, a Taylor series to first order in (u − 1) was used, with u = 1 − iη. Further, an approximation involving η = ω/ω C was made in F . Both steps imply η < 1. The restriction η < 1 means that Eqs. 41 and 57 must be applied to the far-field. Given that a and θ C are fixed by cascade physics and ice density, and that the relevant Askaryan bandwidth for ice is [0.1 − 1] GHz, the parameter most easily varied within η is the observer distance r. Taking ν = 0.5 GHz, n − 1.78, c = 0.3 m GHz, θ = θ C , and a = 5 m, requiring that η = 1 gives r ≥ 0.4 km. Scaling to ν = 0.25 GHz gives r ≥ 0.2 km. According to NuRadioMC [23] (Fig. 13) , the r corresponding to UHE-ν at 10 18 eV ranges from 0.7-3.2 km, and 0.2 km is rare.

The "acceleration argument" invoked by RB in [31] states that if r(t) points to the ICD, r(t) must be constant enough to ensure that ∆r < λ. Using the law of cosines, with two sides being r and r + ∆r, and a third being a, the criteria that (a/r) 2 1 leads to |∆r| ≈ a/n which is O(2) m. When in doubt about usage and event geometry, determining if (a/r) 2 1 is a good check. If the UHE-ν event is a charged-current interaction with an electromagnetic cascade above the LPM energy [20] , a grows faster than ln(E C /E crit ) (Sec. III B). Since future work with this fully analytic treatment will focus on developing probability distributions for parameters like a as part of an energy reconstruction, it will straightforward to study in parallel when the far-field requirement breaks down by comparing to ARZ models in NuRadioMC [23] .

There are at least four main advantages of fully analytic models that accurately model the central Askaryan radiation peak. First, when analytic models are matched to observed data, cascade properties may be derived directly from the waveforms. Second, in large scale simulations required for the design of IceCube-Gen2 radio, evaluating a fully analytic model technically provides a speed advantage when compared to the semi-analytic parameterizations. Third, fully analytic models, combined with RF antenna response, can be embedded in firmare to form a matched filter that enhances detection probability. Fourth, parameters in analytic models may be scaled to produce results that apply to media of different density than ice. This application is useful for understanding potential signals in the Antarctic firn, or the upper layer of snow and ice that is of lower density than the solid ice beneath it.

The ability to fit cascade properties from waveforms will be a useful tool for IceCube-Gen2 radio. Examples of current reconstruction techniques include the forwardfolding method [25] and information field theory (IFT) [26] . In particular, the longitudinal length parameter a leads to a reconstruction of ln(E C ), given knowledge of ∆θ ( Fig. 11 and Equation 71 ). Further, all designs for detector stations in IceCube-Gen2 radio include many distinct RF channels and one phased-array of channels. Treating each waveform like a separate trial would increase cascade parameter precision. Even pulse width measurements without a priori knowledge of a or θ contain information about (θ, a) via the uncertainty principle (see gray contours of Figures 4 and 5) .

Much thought has been devoted to improving the speed of NuRadioMC and other codes that rely on Askaryan emissions models. In one example, AraSim, convolution with pre-simulated cascade profiles was dropped in favor of on-the-fly convolution with the Greisen (EM) and Gaisser-Hillas (hadronic) cascade profiles [22] . While this reduced the number of steps required to produce a simulated data set, it required the code to make a determination of the optimal time-window before convolution. To boost speed in NuRadioMC, a shower library is implemented that pulls realistic cascade profiles from memory and performs the convolutions event-by-event. This strategy provides accuracy regarding randomly forming sub-showers due to the LPM effect, but it still requires a convolution for each event [23] . A fully analytic model provides fast accuracy in the energy range of 10-100 PeV without the need for a shower library.

The most intriguing usage for a fully analytic Askaryan model would be to embed the model as a matched filter in detector firmware. Matched filters are commonly used in radar systems that maximize detected SNR by correlating incoming data with a pre-defined function. Because the cascade properties would not be known in advance, an array of matched filters could be implemented to form a matched filter bank. One example of this approach was the TARA experiment [37] , which was designed to detect low-SNR cosmic ray radar echoes. This is similar to the challenge faced by IceCube-Gen2 radio: pushing the limit of low-SNR RF pulse detection in a remote setting. For example, a matched filter bank could be formed with an array of off-cone field formulas with fixed a-value and varying θ-values, which would then be convolved with the RF channel impulse response (see Section 6 of [27] ).

Finally, a fully analytic model enhances the ability of IceCube-Gen2 radio to identify signals that originate in the firn. At the South Pole, the RF index of refrac-tion begins around 1.35 and does not reach the solid ice value of 1.78 until 150-200 meters [28] . There are at least two signals that could originate in the firn: UHE-ν events that create Askaryan radiation, and UHE cosmic ray cascades partially inside or fully inside the firn. The altitude of the South Pole makes the latter possible. The Askaryan radiation of the firn UHE-ν events could be modeled via appropriate density-scaling of the longitudinal length a and the lateral width 1/( √ 2πρ 0 ). Being lower in density, the firn woud contain cascades of different size and shape relative to solid ice, which would alter the radiation. The Askaryan radiation of UHE cosmic ray cascades within the firn could be similarly modeled, with the caveat that the cosmic ray cascades would also generate geosynchrotron radiation. In the atmosphere, the geosynchrotron radiation dominates. In the firn, the UHE cosmic ray cascade would only be meters in length instead of kilometers, and the Askaryan radiation is expected to dominate.

To begin the evaluation of F , the exponent −i q · x must be expanded. Let nk → k, and assume that k now refers to the scalar wavevector in the medium.

From the coordinate system, ρ = r sin(θ), so

Then we have for F : (A4) Two physical parameters are necessary to compare wavelength to lateral ICD width. First, let γ be the projected scalar wavevector, and σ be the ratio between γ and the ICD exponential scale factor:

The integral becomes

Rearranging to perform the φ -integral first:

The cylindrical symmetry of f (x ) leaves one free to choose where φ = 0. Thus, the azimuthal angle φ → φ − π may be rotated, changing − cos(φ ) to cos(φ ):

The φ -integral is now proportional to one representation of the 0-th order Bessel function of the first kind: 

According to Equation A17, the form factor is a lowpass filter with the cutoff-frequency ω CF introduced in Section II: σ = ω/ω CF .

(1 + (ω/ω CF ) 2 ) 3/2 ≈ 1 + 3 2

ω ω CF 2 (B7)

Using the approximations introduces simple poles into the complex formula for the frequency-dependent electric field. Inserting the approximations in the denominator of Eq. B5, we have r E(ω, θ C ) = (−iω)E 0 sin(θ C )e iωR/c 1 − i 2 ω/ω C 1 + 3 2 (ω/ω CF ) 2 (B8)

The denominator can be rearranged by factoring the ω coefficients, and defining ω 0 = 2 3 ω CF .

r E(ω, θ C ) = 2iω C ω 2 0 (−iω)E 0 sin(θ C )e iωr/c (2iω C + ω) (ω + iω 0 )(ω − iω 0 ) (B9)

LetÊ 0 = E 0 sin(θ C ), and let the retarded time be t r = t − r/c. Taking the inverse Fourier transform, using the same sign convention as RB [31] (f (t) = (2π) −1 ∞ −∞ F (ω)e −iωt dω), converts the field to the timedomain:

1. If t r > 0: Consider the contour comprised of the real axis and the clockwise-oriented negative infinite semi-circle. On the contour, the exponential phase factor in Eq. B10 goes as exp(−iωt r ) = exp(−i(R cos φ + iR sin φ)t r ) (B11)

For the semi-circle, φ ∈ [π, 2π], so sin φ < 0 and t r > 0. Exponential decay occurs and the integrand vanishes on the semi-circle for |ω| = R → ∞.

2. If t r < 0: Consider the contour comprised of the real axis and the counter-clockwise-oriented positive infinite semi-circle. On the contour, the exponential phase factor in Eq. B10 goes again as exp(−iωt r ) = exp(−i(R cos φ + iR sin φ)t r ) (B12)

For the semi-circle, φ ∈ [0, π], so sin φ > 0 and t r < 0. Exponential decay occurs and the integrand vanishes on the semi-circle for |ω| = R → ∞.

Using cases 1 and 2, Equation B10 can be solved using the Cauchy integral formula. Beginning with t r > 0, two poles are enclosed in the semi-circle: one that originated from the coherence cutoff frequency, and the other that originated from the form factor. The Cauchy integral formula yields

Define the ratio of the cutoff frequencies: = ω 0 /ω C . After evaluating the time derivatives, Equation B13 becomes rE(t, θ C ) =Ê 0 ω 2 0 e −2ωCtr (1 − 2 )(1 + 2 ) − e −ω0tr (2)(1 − 2 ) (B14) Expanding to linear order in , assuming < 1, and recalling that ω 2 0 = 2 3 ω 2 CF :

(B15) Turning to the case of t r < 0, consider integrating Eq. B10 along the contour comprised of the real axis and the counter-clockwise-oriented positive infinite semicircle. The contour encloses one pole, and the exponent ensures convergence: rE(t, θ C ) = (2πi)Ê 0 (π) −1 iω C ω 2 0 d dt r e ω0tr (2iω C + iω 0 ) (2iω 0 ) (B16) After evaluating the derivative, the expression simplifies with = ω 0 /ω C :

Finally, using the same first-order approximation in as the t r > 0 case:

Collecting the t r > 0 and t r < 0 results together:

e ω0tr t r < 0 2e −2ωCtr − 1 + 1 2 e −ω0tr t r > 0 (B19)

Using Tabs. II-IV, Eq. B2 reduces to

Expanding to first-order with respect to u near (u = 1) gives E(u, x) = E(x, 1) + (u − 1)Ė(x, 1) + O(u − 1) 2 (C2)

The first term is f g(1−h) evaluated at u = 1: exp(−y) ( Table IV) . The second term requires the first derivative of E(u, x) with respect to u, evaluated at u = 1.

E(u, x) = fġ +ḟ g − (f gḣ + fġh +ḟ gh) (C3) E(1, x) = fġ +ḟ g − (f gḣ + fġh +ḟ gh) | u=1 (C4)

The first-derivatives of f , g, and h, evaluated at u = 1, are given in Tab. IV. Because h(x, 1) = 0, terms proportional to h will vanish. The result iṡ

Inserting Eq. C5 into Eq. C2, E(u, x) = e −y 1 + 1 2 (u − 1) (2y + 2q − 1) (C6)

Using the definition of u (Table II) , the result may be written E(u, x) = e −y 1 − 1 2 jη (2y + 2q − 1) (C7)

Proceding with the inverse Fourier transform of thê θ-component:

Let η = ω/ω C , y = pω 2 (Table II) . Inserting the Taylor series for E, the form factor F , and ψ = −i exp(ikr) sin θ (Sec. II), and following the same steps as the on-cone case produces

Unlike the on-cone case, Equation C9 cannot be in-tegrated with infinite semi-circle contours, because the exponential term diverges along the imaginary axis far from the origin. Let I 0 represent the constant term with respect to ω in the numerator:

Further, let I 1 and I 3 represent the linear and cubic terms, respectively. Completing the square in the exponent of I 0 , with ω 1 = t r /(2p), yields

may be re-cast as the line-broadening function, H (DLMF 7.19, [34] ) common to spectroscopy applications:

Assume that ω > ω 1 . This approximating step will be called the symmetric approximation.

The result for I 0 involves the complementary error function (DLMF 7.7.1, [34] ):

2p πω −1 0 e pω 2 0 erfc( √ pω 0 ) (C14)

The integrals I 1 and I 3 are zero by symmetry, with odd integrands over (−∞, ∞). Inserting the result for I 0 into Eq. C9 and evaluating the derivative finishes the problem (see Sec. V).

Radioglaciology

Release 1.1.1 of 2021-03-15

Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables

We would like to thank our families for their support throughout the COVID-19 pandemic. We could not have completed this work without their help. We would also like to thank our colleagues for helpful discussions regarding analysis techniques. In particular, we want to thank Profs. Steve Barwick, Dave Besson, and Christian Glaser for useful discussions. Finally, we would like to thank the Whittier College Fellowships Committee, and specifically the Fletcher-Jones Fellowship Program for providing financial support for this work. This work was partially funded by the Fletcher-Jones Summer Fellowship of 2020, Whittier College Fellowships program.

The original equations for theθ-component of E are:The complete field from the original RB model [31] , including the form factor F , ψ = −i exp(ikr) sin θ, and E isLet Eq. 8 for the form factor, with σ = ω/ω CF and η = ω/ω CF , and letting E 0 be proportional to cascade energy E C :Suppose ω < ω C , and ω < ω CF , such that the following approximations of the factors in the denominator are valid: