Hypothesis Tests on Rayleigh Wave Radiation Pattern Shapes: A Theoretical Assessment of Idealized Source Screening

Shallow seismic sources excite Rayleigh wave ground motion with azimuthally dependent radiation patterns. We place binary hypothesis tests on theoretical models of such radiation patterns to screen cylindrically symmetric sources (like explosions) from non-symmetric sources (like non-vertical dip-slip, or non-VDS faults). These models for data include sources with several unknown parameters, contaminated by Gaussian noise and embedded in a layered half-space. The generalized maximum likelihood ratio tests that we derive from these data models produce screening statistics and decision rules that depend on measured, noisy ground motion at discrete sensor locations. We explicitly quantify how the screening power of these statistics increase with the size of any dip-slip and strike-slip components of the source, relative to noise (faulting signal strength), and how they vary with network geometry. As applications of our theory, we apply these tests to (1) find optimal sensor locations that maximize the probability of screening non-circular radiation patterns, and (2) invert for the largest non-VDS faulting signal that could be mistakenly attributed to an explosion with damage, at a particular attribution probability. Lastly, we quantify how certain errors that are sourced by opening cracks increase screening rate errors. While such theoretical solutions are ideal and require future validation, they remain important in underground explosion monitoring scenarios because they provide fundamental physical limits on the discrimination power of tests that screen explosive from non-VDS faulting sources.


Introduction
Seismic sources produce Rayleigh wave radiation patterns that indicate focal mechanism asymmetry (Fig. 1). Teleseismic and regional records of large underground nuclear explosions (∼ 10 1 − 10 3 kT) emplaced in pre-stressed geologies, for example, have historically revealed that tectonic release of strain can superimpose with the isotropic explosion source to produce Love and Rayleigh waves with non-circular radiation patterns (Ekström and Richards, 1994;Toksöz and Kehrer, 1972). These observations contrast with simple models of cylindrically-symmetric, shallow sources hosted in vertically stratified geologies that predict azimuthally constant Rayleigh wave amplitudes and absent Love waves. Instead, such observations are partially consistent with models of faulting motion presented by tectonic release, which produces asymmetric deformation about a vertical axis that excites a Rayleigh-wave field with a non-circular radiation pattern (Richards and Ekström, 1995). Observers risk misidentifying large, shallow explosions that accompany such tectonic release as shallow earthquakes when the imprint of such a "faulting signal" on the radiation pattern shape is significantly non-circular or has a low signal-to-noise ratio (SNR). Smaller yield explosions (∼ 10 1 − 10 3 kg) that do not trigger tectonic release may still accompany fracturing events or mechanical sliding on rock joints that output shear energy as bulk damage (Steedman et al., 2016). These bulk effects do not appear to significantly distort surface wave radiation patterns when they are symmetric about a vertical axis. In particular, local distance records from a series of conventional explosives conducted through the Source Physics Experiment (SPE) (Snelson et al., 2013) demonstrate that low yield, shallowly buried sources that cause damage are vertical-axis symmetric, and thereby trigger Rayleigh wave radiation patterns that are circular within a factor of 1.25 (Larmat et al., 2017). These explosions, which were emplaced in the granites of Climax stock on the Nevada National Security Site (NNSS), also demonstrate that body waves radiated by the same source in a similar frequency band show apparent azimuthally asymmetry that is likely sourced by refraction effects (Darrh et al., 2019). Experiments with low-yield explosives emplaced in boreholes within New Hampshire granite demonstrate more complex effects. Namely, these latter experiments reveal that observationally confirmed radial fracturing in hard rock accompanies asymmetric, short period (∼6Hz), guided Rayleigh waves (Rg). Body waves sourced by these same New Hampshire explosions can exhibit even greater azimuthal dependence (Stroujkova, 2018). Both the SPE and New Hampshire experiments cumulatively indicate two substantive effects of the seismic source on radiation pattern shape can be observed from local distances. First, vertical axis-symmetric damage does not significantly  Figure 1: Surface sensors (triangles) record Rayleigh waves from a shallow explosion that is buried to a depth h in a horizontally-stratified half-space that has moment M I , which creates damage with moment M CLVD that each superimpose with faulting of moment M 0 . The damage occurs when constructive interference between the incident and reflected waves create a conicallyshaped zone of failure. (a): An explosion source capable of producing damage that is symmetric about a vertical axis, but that includes no faulting (no tectonic release) radiates a circular radiation pattern R EX described at top. (b): A source that includes non-VDS faulting near the explosion point radiates a non-circular radiation pattern R EQ described at top. Each radiation pattern shape associates with a deterministic scalar Λ (bottom) that quantifies the power of the hypothesis test between data that record circular versus non-circular radiation patterns (at top). Subsection 4.1 defines each symbol. distort Rayleigh wave radiation pattern shapes from circularity, whereas a structure at depth can distort body wave patterns. Second, significant radial fracturing and shearing motion at the shallow explosion source measurably distorts both Rayleigh wave and body wave radiation patterns. These observations indicate that azimuthal variability of Rayleigh waves sourced by small yield explosions indicate faulting or shearing motion at the source, whereas body wave asymmetry appears non-uniquely attributed to source and path.
A forthcoming test within the SPE series will provide a unique opportunity to study asymmetry in radiation patterns that are sourced by explosions and superimposed with historic, tectonic records of fault slip. This effort will detonate a conventional explosive in tectonic regions and monitor its seismic signatures for comparison against earlier earthquake signatures that output radiation that shares wavefront paths (Walter et al., 2012). If Rayleigh waves output by this explosion remains circular and mismatches historical records, theory and the aforementioned observations imply that the resultant radiation pattern shape provides a hopeful discriminant. It remains unclear, however, "how much" non-circularity of the Rayleigh wave radiation pattern is sufficient to indicate that a seismic source is dominated by shearing or fault motion, and is therefore more earthquake-like than explosion-like. A discriminant that tests radiation pattern distortion from local to teleseismic ranges can therefore benefit both the source physics community and support nuclear test-ban treaty surveillance, if justified on both theoretical and experimental grounds.
This work determines if a discriminant for shallow source types that tests Rayleigh wave radiation pattern shapes is justified on theoretical grounds. To confront this challenge, we apply binary hypothesis tests to competing data models of noisy radiation patterns and assume several source parameters are unknown. Our models also use planar geologies and therefore isotropic propagation of Rayleigh waves; we do not consider Love waves. Under this limited scope, we derive screening curves from the hypothesis tests that bound an observers's capability to discriminate non-vertical dip-slip (non-VDS) shallow faulting sources from symmetric sources (like shallow explosions that create damage that is symmetric about a vertical axis). The screening power of these tests is completely defined by the product of a faulting-signal term, and a term that depends on deployment parameters. Using this test, we illustrate several applications of our theory that include finding deployment locations for sensors to supplement an existing network and that maximize our source-type screening probability. We conclude that a Rayleigh wave radiation pattern discriminant shows a high probability Pr D of success (Pr D > 0.9) when a sufficient number of sensors (> 12) with a limited azimuthal gap (< 90 • ) records radiation from sources with moderate strike-slip faulting signals (SNR > 20). This probability diminishes with more ambiguous focal mechanisms that resemble certain historical events recorded near the Lop Nor nuclear test site in China.
Our results suggest that a screening statistic that tests Rayleigh wave radiation pattern shapes for faulting signals requires very dense sensor coverage to achieve good efficacy, and is probably impractical in most passive monitoring scenarios. We lastly consider how radiation pattern distortion from opening cracks and wavefield homogenization can inflate screening errors.
We emphasize that our treatment is purely theoretical and idealized. Our models require that future efforts provide any evidence that this approach has a practical utility as a screening technique in non-ideal settings. In particular, focusing effects (Lin et al., 2012;Selby and Woodhouse, 2000), tectonic pre-stress (Ekström and Richards, 1994;Toksöz et al., 1971), surface roughness (Steg andKlemens, 1970), path-dependent attenuation (MacBeth andBurton, 1987), transmission losses through vertical interfaces (Napoli and Russell, 2018), and scattering from topography that is comparable in relief to Rayleigh wavelengths (Ichinose et al., 2017;Snieder, 1986) often dominate the azimuthal variability of observed surface wave radiation. Stratified half-space models that support idealized Rayleigh wave propagation are likely limited to very local seismic measurements of SPE shots, cyroseismic studies interior to crevasse-free ice sheets and glaciers (Carmichael et al., 2015;Hudson et al., 2020;Lindner et al., 2019;Lough et al., 2015), and to engineered structures (Lu et al., 2007). We assert that our study remains useful in geophysics, however, because it helps define the physical limits on an observer's capability to screen Rayleigh wave sources, using only their radiation pattern shapes.
Such limits support previous efforts that exploit Rayleigh wave radiation pattern geometry to better characterize underground explosions (Toksöz et al., 1971;Yacoub, 1981), and supports present computational and observational work to understand earthquake radiation (Rösler and van der Lee, 2020).

Reference Sources and Historical Events
Our assessment will compare cylindrically symmetric sources that include explosions with damage along the vertical axis against a set of four distinct non-VDS, tectonic earthquake sources. One such faulting source is the strike-slip earthquake. This particular double-couple source is maximally dissimilar from any isotropic source, as plotted Lop Nor (2003) Figure 2: Three focal mechanism solutions for shallow earthquakes that locate near underground test sites and the noisy, non-circular component of the radiation pattern for each source. (a): The M W = 3.7, 1993-05-30 event from the Rock Valley Sequence, termed "Event 1" (depth ≈ 2 km) (Smith et al., 1993). (b): The M L = 3.3, 1997-08-16 event from the Kara Sea (0 km < depth < 23 km) that located tens of km from Novaya Zemlya (Bowers, 2002;Hartse, 1998;Schweitzer and Kennett, 2007). c: The m b = 4.8, 2003-03-13 tectonic event located near the Lop Nor test site (5 km < depth < 7 km) (Selby et al., 2005). In each pattern, Gaussian noise superimposes with a deterministic radiation pattern solution to mimic a finely sampled observation with an overall SNR of 20. We note that the noise on the non-circular component of the radiation pattern increases from left to right, as the sources produce Rayleigh waves less efficiently. We further assume that the isotropic component of the source is zero in each case, and within the bound of any uncertainties in the focal mechanism estimate. Here, M W is the moment magnitude, M L is local or Richter magnitude, and m b is body wave magnitude.
on a Hudson diagram (Hudson et al., 1989) or lune (Tape and Tape, 2019

Shallow Source Rayleigh Waves
We consider Rayleigh wave motion present in a cylindrically symmetric, horizontally layered half-space triggered by a buried seismic source (Fig. 1). This source may be described by a point force or a symmetric moment 3 × 3 moment tensor M that summarizes displacement discontinuities as sets of force couples. Such a moment tensor source excites the frequency-domain, Rayleigh wave displacement u = [u x , 0, u z ] T that solves Newton's second law of motion in this stratified half-space, subject to traction-free surface boundary conditions. The scalar product of M with the gradient of the (Hermitian) Rayleigh wave Green's tensor G RAY (Aki and Richards, 2002, eq. 3.23) compactly represents such a general displacement: When the source is emplaced in half-space with non-trivial stratification (not pure halfspace), the displacement solution of eq. 1 is a linear combination of eigenvectors r that sum over distinct wavenumber modes (Aki and Richards, 2002, eq. 7.150-7.151). The dispersion relationship k = k (ω) for this Rayleigh wave field defines these modes and couples the radiation pattern R of u to the range-and depth-dependent part of the displacement solution g. If the moment tensor source is shallow enough that the freesurface traction free boundary conditions σ ·ê z = 0 also apply at the source's depth, then the radiation pattern effectively decouples from g (Aki and Richards, 2002, pg. 328); here, σ is the elastic stress tensor,ê z is a unit vector pointing vertically upward, and 0 is a three element vector of zeros. More formally, this condition requires that the source is buried at a depth that is small compared to the dominant wavelength of the seismic radiation triggered by the source-time function. These cumulative conditions mean the Rayleigh wave field has the simple form and is a product of R and g: Vector g in eq. 2 depends on source depth and observation range, but not φ; termR is the mean radiation pattern; R m (m = 1, 2) are non-circular radiation pattern coefficients that all have units of moment and depend on moment tensor M ; and φ is the azimuthal angle (see Aki and Richards, 2002, Fig. 4.20). The general displacement in eq. 1 and the traction free boundary conditions relate the mean radiation pattern and azimuthally dependent coefficients to linear functions of the moment tensor elements: where M ij is the force couple aligned with the Cartesian direction i separated by direction j. We note that Rayleigh wave radiation patterns do not resolve terms M xz and M yz that indicate vertical shearing motions directed parallel to Earth's surface.

Colocated Faulting, Explosion and Damage Sources
The most general symmetric moment tensor for a single point source is a superposition of an explosion with isotropic moment M I that effectively colocates with a non-VDS faulting source of moment M 0 and a damage source, which is equivalent to a compen- drive contraction along that axis. We assume that the fault is parameterized by strike (φ s ), rake (λ) and dip (δ) angles. The radiation pattern coefficients from eq. 2 are then (algebra omitted):R = 2β 2 where we have used a moment decomposition (Patton and Taylor, 2011, eq. 15). In eq. 4, scalar SS quantifies the strength of the strike-slip faulting component, whereas scalar DS quantifies the strength of the dip-slip component. These terms are (Ekström and Richards, 1994, eq. 21 and eq. 22): where the equivalent expression for SS elsewhere (Patton and Taylor, 2011, eq. 16) is missing a factor of two in the dip angle. We express the full Rayleigh wave displacement by combining eq. 2, eq. 4, and eq. 5. The superposition of the explosive, damage and non-VDS faulting sources thereby creates a Rayleigh wave displacement field that is: where ψ = φ − φ s is the difference between the azimuthal and strike angles; because we use ψ hereon, we refer to it as the "azimuthal" angle, despite abuse in terminology.
Regardless, the function that scales g is the radiation pattern R, and includes the same functional form as the scalar factor that appears in eq. 2: Eq. 4 with eq. 7 demonstrates thatR is azimuthally constant, and that azimuthal variability arises from non-VDS faults where δ = ±π/2, λ = ±π/2. Eq. 7 is also the starting point for estimating the optimal receiver sampling of the Rayleigh wave radiation pattern of shallow sources (Section 4.3).

Statistical Hypothesis Tests
To determine if observed Rayleigh wave radiation indicates an isotropic or mixed non-VDS source, we evaluate a hypothesis test that compares two distinct models for deterministic signals R. This test forms a generalized likelihood ratio statistic from N observations of R at different azimuthal locations and compares them to a threshold to measure evidence that the data record a non-circular radiation pattern. Specifically, our null hypothesis is that an azimuthally distributed set of sensors samples a circular radiation pattern with zero non-VDS faulting moment (M 0 = 0). Our alternative hypothesis states that these same sensors sample a non-circular radiation pattern. The competing hypotheses are (without noise): where δ, λ, and M 0 parameterize DS and SS, and are unknown. We now add a random variability models to eq. 8 in two respects: first, we superimpose Gaussian measurement noise to R, and second, we model R itself to be a sum of the deterministic component and a zero mean random process. To include measurement noise, we note that observations often show that moment estimates are log-normally distributed if such moment is measured from station magnitudes (Godano and Pingue, 2000;Kagan, 2002). Moment estimates that fit a horizontal line to the low-frequency component of waveform's displacement spectrum reveal Gaussian residuals, because normally distributed errors dominate seismic waveform records (Šílenỳ et al., 1996;Walter et al., 2009). We therefore assume waveform noise is the dominant source of random error in scalar moment estimates and model a single observation at fixed azimuth as a Gaussian random variable.
To modify the signal portion of the radiation pattern, we add a zero-mean Gaussian process with covariance σ 2 R to the deterministic, non-zero mean of R. We indicate that R is a random variable with mean µ and variance σ 2 with notation R ∼ N µ, σ 2 , where ∼ indicates "is distributed as". Eq. 8 then tests between two competing probability distributions for R (with noise): Eq. 9 applies to one observation ψ. For many azimuthally distributed observations, our hypothesis test includes a vectorized set of radiation pattern samples R where the k th sample is R k : The competing hypotheses (eq. 9) form a general Gaussian detection problem, applied to a multi-sample radiation pattern (Kay, 1998, section 5.6). To explicitly define each radiation pattern term, we concatenate the unknown statistical parameters into a single vector. We then express the mean µ of R under H 1 as a linear system that equates with eq. 4 so that component m (row m) of µ is: where: and where H(m, :) in eq. 11 defines row m of matrix H. The first parameter ∆R weights the circular portion of the radiation field whereas the other parameters weight the noncircular portion. The covariance between the observations R i and R j at two azimuths is zero (eq. A.5). This independence in observations implies that the probability density function (PDF) for vector R under H 1 has mean µ = Hθ and covariance (σ 2 R + σ 2 M )I. The PDF under H 1 is therefore: where I is the identity matrix. Under H 0 , DS = SS = σ 2 R = 0. The mean µ of R under H 0 is also linear system, but θ is now constrained to make R circular: The PDF for vector R with covariance σ 2 M I under H 0 is then: These competing PDFs provide the basis for a binary hypothesis test between circular versus noncircular radiation patterns. The hypothesis test in eq. 8 that uses such vectorized data is: To exploit these tests, we consider two cases. This first case (Case I) considers the variability σ R in R from random processes as known, but the stochastic variability σ M as unknown. The second case (Case II) considers σ M as approximately known, but σ R as unknown, though in much less detail than Case I. Parameter θ is unknown in both cases. Case I and Case II both include the scenario that σ R is zero as a special example.
4.1 Case I: σ R known, θ and σ M unknown We first assume a network of N ≥ 8 sensors records a Rayleigh wave radiation pattern with known process variance (σ 2 R in eq. 16), unknown faulting parameters (θ in eq. (Case I) Figure 3: Probability density functions that describe the source screening statistic L θ (R) (horizontal axis) in eq. 18. The statistic describes N = 10 receivers randomly distributed over a circle centered at the hypocenter of a shallow seismic source (orange star) that sample a source radiation pattern R. The leftmost thick, gray curve indicates a circular radiation pattern (Λ = 0) and rightmost curves indicate noncircular radiation patterns (Λ > 0). The darkest shaded area under the null PDF curve (Λ = 0) indicates a Pr F A = 5 · 10 −3 false screening probability. The lighter shaded area below the alternative PDF (Λ = 10) measures the screening probability Pr D = 0.81. The threshold is η ≈ 5 (eq. 23). 11), and unknown stochastic variance (σ 2 M in eq. 16). This case idealizes explosion monitoring scenarios in which an underground testing complex colocates with a shallow fault of unknown geometry. This case also includes the scenario that all variability is stochastic in origin (σ 2 R = 0) as an example. Under this general model, a sensor network then records a low magnitude seismic event that originates near a test site, and an observer must estimate the probability that this event is explosive or tectonically sourced by a non-VDS fault (hypothesis H 0 or H 1 ).
To determine which hypothesis Rayleigh wave records likely describe, we use a generalized log-likelihood ratio test (log-GLRT). This test compares the logarithmic ratio of the alternative hypothesis PDF to the null hypothesis PDF, where each density is evaluated at maximum likelihood estimates (MLE) of its unknown parameters to form a test statistic. The log-GLRT is equivalent to a conditional decision rule on this scalar test We compute the MLEs of the parameters in Appendix B using projector matrices and reduce the scalar statistic into a ratio of two scaled, statistically independent terms. This provides a test statistic L θ (R) equivalent to eq. 17 for screening circular from non-circular radiation patterns: The rank-1 projector matrix P X in the numerator of eq. 18 is: The rank-3 projector matrix P H = H H T H −1 H T projects vectors onto the subspace spanned by H, span{H}. Projector matrix P ⊥ H = I -P H in the denominator of eq. 18 then projects vectors onto the space orthogonal to span{H}. We thereby interpret the screening statistic L θ (R) in eq. 18 as a scaled ratio of the radiation pattern energy due to non-VDS faulting (note structure of A), divided by the residual radiation pattern energy not attributable to the Rayleigh wave model (the system H). Probability theory predicts that this quotient has a noncentral-F distribution with one, and N − 3 degrees of freedom. The density for L θ (R) is then shaped by a scalar noncentrality parameter Λ if R k includes Gaussian noise (1 ≤ k ≤ N ) and the radiation pattern is non-circular. We write its distributional dependence as L θ (R) ∼ F 1,N −3 (L, Λ > 0) and the cumulative, noncentral F distribution for L θ (R) as: These distributions have a well-defined variance when their second degree of freedom exceeds four, so that the screening statistic L θ (R) is only well-defined when N ≥ 8. This sampling implies that at least two sensors sample each quadrant of the Rayleigh radiation pattern (given full azimuthal coverage by sensors), which mitigates spatial aliasing of four-lobed radiation patterns. This statistic is identical in form to the subspace detector (Scharf and Friedlander, 1994) that is often applied in seismic monitoring, and has completely analogous interpretations for waveform detection. Consequently, many standard detection theory texts document the properties and analytical form for F 1,N −3 (L, Λ), and computational packages like MATLAB include routines to calculate random variables and parameters of the noncentral-F distribution.
Crucially, this theory shows that the scalar Λ in eq. 20 completely quantifies the screening capability of L θ (R) to test between H 0 and H 1 at fixed N ; eq. B.5 shows its general analytical form. Here, we use those results to write Λ as the product of one factor that depends on deployment azimuth ψ and a second factor that depends on source and noise parameters: This factorization shows that Λ is the product of a scalar term A H T H −1 A T −1 that stores the N sensor deployment locations, and a term (DS + SS) 2 /(σ 2 M + σ 2 R ) that quantifies the non-VDS faulting signal energy, relative to noise and random process variance. It is analogous to the signal-to-noise plus interference ratio (SNIR) of waveform processing. Explosive or VDS sources that trigger no Rayleigh waves include parameters σ 2 R = DS = SS = 0, and Λ = 0 (as expected). Strike-slip faults maximize Λ among other faulting sources, given equal sensor placement, because (DS + SS) 2 / σ 2 M + σ 2 R ≤ M 2 0 / σ 2 M + σ 2 R (see eq. 5). Importantly, eq. 18 and eq. 21 demonstrate that discriminating between a circular and non-circular Rayleigh wave radiation pattern under the Case I assumptions is entirely equivalent to a signal detection problem. We further note the random process variance σ 2 R acts to simply inflate the effective noise variance and reduce any observed, non-VDS faulting signal SNR under H 1 . This asymmetry in variance between the null and alternative hypotheses means that a high SNR faulting signal with high random signal variance is indistinguishable from a lower SNR faulting signal with zero random signal variance. We will assume σ 2 R = 0 hereon, so that Λ and L θ (R) are each simplified from their more general forms. This simplification preserves our physical insight into the radiation pattern testing problem without adding cumbersome notation. The faulting SNIR then reduces to just SNR.
Under these simplifying (but still rather general) assumptions, the probability Pr D of detecting any general faulting signal, and then deciding that R records a non-circular radiation pattern is the probability Pr D of measuring a value of L θ (R) that exceeds a threshold η (eq. 18): We select a threshold η for deciding a radiation pattern is non-circular using the Neyman-Pearson criteria. This criteria uses a constant false-attribution probability Pr F A to invert for η: where we set Pr F A to an acceptable value (e.g., 10 −3 ) and F −1 1,N −3 (•) is the inverse of the central F distribution with one and N − 3 degrees of freedom. The screening probability Pr D in eq. 22 monotonically increases with Λ at fixed N and η. The locus of points (Λ, Pr D ) thereby defines performance curves that measure the screening capability of the decision rule in eq. 18 versus SNR of the faulting signal at fixed Pr F A (Fig. 3).
Deployment constraints and tectonic settings can each limit the permissible parameter space of Λ to a subset S (Λ). These limits include the azimuthal deployment range of sensor numbers and locations that quantify the leftmost factor in eq. 21, and variability in the source focal mechanism, which quantifies the rightmost factor in eq. 21. In such cases, we estimate the expected screening capabilityPr D of eq. 18 as the average of the detection rate Pr D over such parametric constraints. We use eq. 22 to write this average as:P in which |S (Λ)| is the size of the set (e.g., the number of realizations in a Monte Carlo simulation). The sum in eq. 24 equates to an integration when the subset of permissible values of Λ is a continuum (e.g., |ψ| ≤ 3π/4) and gives a marginal distribution forPr D (e.g.,Pr D = 4/ (3π) 24) that each index distinct deployment and source parameters. Fig. 4 shows the probability that eq. 24 will screen strike-slip earthquakes from explosions with N = 12 sensors deployed over three distinct azimuthal ranges against faulting signal (pictured at right). Each of these three curves show the dependence of screening probability on the faulting signal SNR (DS + SS) 2 /σ 2 M , averaged over 100 individual screening curves. Scalar Λ uniformly samples receiver deployments confined to these azimuthal ranges where source parameters are δ = π 2 , λ = 0, (DS + SS) 2 = M 0 , and faulting SNR = M 2 0 /σ 2 M . The three curves demonstrate that screening probability markedly decreases as azimuthal gap increases with a fixed number of stations, even as station density over deployment azimuth (sensor number per arc length) increases. Shading about the top blue curve associated with full azimuthal sampling marks the one-standard deviation (±σ) variability. Fig. 5 illustrates that earthquake signals triggered by other focal mechanisms are more difficult to screen from explosions relative to strike-slip events (δ = ±π/2, λ = 0, π).
Specifically, eq. 18 illustrates predicted curves that quantify the capability of statistic L θ (R) to screen distinct earthquakes from explosions using the same deployment strategy with N = 12 sensors. In this case, each curve is an average of 100 screening curves that randomly sample a circular instrument deployment around three different sources (focal mechanisms at right). These curves demonstrate that the radiation pattern test screens non-strike-slip earthquakes from explosions less effectively than strike-slip earthquakes, even with identical station coverage. Seismic events with a Rayleigh wave radiation pattern like that of the 2003 Lop Nor tectonic event that is sampled by an unrestricted deployment of N = 12 sensors is particularly difficult to screen from an explosion. Shading about the curve that associates with full azimuthal sampling marks the ±σ variability in screening sources like Event 1 from the Rock Valley earthquake. Figure 4: Averaged screening curves (eq. 24), plotted against faulting signal SNR (horizontal axis; eq. 21). Each curve presents an estimate of 100%× the probability Pr D (vertical axis) that the decision rule in eq. 18 will correctly screen a strike-slip earthquake (beach ball, top left) from an explosion. The test statistic L θ (R) (eq. 18, left-hand side) for each curve associates with N = 12 sensor azimuths that sample the Rayleigh wavefield R under distinct deployment constraints (right); all thresholds η correspond to Pr F A = 10 −3 (eq. 23). The topmost curve (blue) presents an average of |S (Λ)| = 100 individual decision rule curves that each measure screening performance when N = 12 sensors uniformly sample random azimuths along the right, topmost dashed circle that surrounds a source (orange star). The shaded region indicates ± one standard deviation estimate (±σ) from the mean curve. The second curve (red) presents a similar average of 100 curves that associate with sensors deployed with uniformly random azimuth along the red dashed arc at right surrounding the same source; in this case, the azimuthal range is restricted to 2π/3 radians. The third curve (green) shows a similar average, but for sensor azimuths restricted to π/2 radians (green dashed arc). The filled circles superimposed on the dashed deployment arcs at right show particular, random sensor locations in each case. Some circular sensor markers overlap.  Fig. 4, but curves now present 100%× the probability Pr D (vertical axis) that eq. 18 will screen earthquakes with three distinct focal mechanisms from an explosion, using N = 12 sensors, all plotted against faulting signal SNR (horizontal axis; eq. 21). Each curve is an average of 100 screening curves (eq. 24) that associate with sensor deployments that uniformly sample all azimuths along the dashed circle centered at the hypocenter of a shallow source (orange star); the filled circles mark a particular, random deployment. The topmost curve (blue) presents such an average for faulting source like that of "Event 1" from the Rock Valley earthquake; the second curve (red) presents a similar average for faulting source like that of the Kara Sea event; and the bottom curve presents a similar average for a faulting source like that of the 2003 Lop Nor event. The shaded region indicates ± one standard deviation from the top curve mean. All curves correspond to Pr F A = 10 −3 . Figure 6: Similar to Fig. 4, but curves now estimate 100%× the probability that eq. 18 (vertical axis) will screen earthquakes whose focal mechanisms match that of the 2003 Lop Nor event (beachball, upper left) from an explosion, using three distinct sensor deployments without azimuthal constraint. As in Fig. 4 and Fig. 5, each curve presents an average of 100 curves (eq. 24) that associate with sensors deployed with uniformly random azimuth along the dashed circles (circular markers, right) surrounding the source (orange star). All curves correspond to a Pr F A = 10 −3 false attribution probability and plot against faulting signal SNR (horizontal axis; eq. 21). Filled sensor circles mark a particular, random deployment.
An increase in sensor deployments (N > 12) correspondingly increases screening probability between non-VDS faulting and explosion sources. Fig. 6 shows three distinct averages of 100 screening curves that each associate with three distinct sensor deployments (N = 8, 16, 24). Each deployment randomly samples a circle centered about a shallow source like that of the 2003 Lop Nor tectonic event. In this case, a screening statistic that uses N = 24 sensors can screen non-VDS faulting from explosion sources at roughly the same mean rate as a statistic that uses 12 receivers to screen sources like Event 1 from the Rock Valley earthquake sequence (Fig. 5). While a 24 sensor deployment screens shallow sources like that show in Fig. 6 with relatively low variance, the test statistic obviously requires twice the number of observations.

A Note on Screening Curve Interpretation
Our interpretation of screening curves that compare probabilities likePr D against faulting signal SNIR require comment (see Fig. 4, Fig. 5 and Fig. 6). The binary decision rules necessarily decide between a continuum of source types that are not identically explosions or faults. For example, the decision rule in eq. 18 will screen data that records Rayleigh motion sourced by an explosion that superimposes with non-zero contributions from a non-VDS faulting signal with a probability that much greater than the falseattribution probability Pr F A = 10 −3 . A proper false attribution probability constraint that considers an acceptable amount of faulting signal with an explosion would more appropriately marginalize out some faulting effects with a prior on Λ from this continuum (Berger and Delampady, 1987). We address such source-type errors and ambiguities in Section 5.4 of our Discussion. The binary hypothesis testing theory that we use here is standard, however, and we proceed with the caveat that neither the H 0 nor H 1 in our models fully capture the continuum of source type hypotheses.
We next explore if more optimal sensor placement can improve screening rates of target sources that are particularly challenging to screen from cylindrically symmetric sources.

Case I application: optimal deployment of auxiliary receivers
We suppose an existing N ≥ 8 receiver network with limited azimuthal coverage requires an auxiliary receiver to maximize its capability to screen shallow non-VDS sources from explosions, through application of L θ (R) as a discriminant. To develop an objective function for this discriminant, we add a row [1, cos(2ψ) − c, sin(2ψ)] to H and make a new matrix H +1 . We then maximize Pr D over additional sensor deployment azimuths ψ +1 that we estimate asψ +1 . Such an optimal auxiliary sensor location estimate supplements the existing network to form an updated N + 1 network that maximizes the probability of detecting a non-VDS signal in the radiation pattern data R. Specifically, this new observation point updates the discriminant in eq. 18, while eq. 23 updates the threshold η to maintain the same false attribution probability Pr F A . The optimal solutionψ +1 depends on both η and Λ, and maximizes Pr D where the total differential dPr D peaks:ψ The decision threshold decreases by an increment ∆η: whereas the faulting signal changes by increment ∆Λ: The argmax argument on the right-hand-side eq. 25 is a directional derivative. It is maximal at parameter increments [∆η, ∆Λ]  since H T H −1 is positive definite. This collective system structure implies that the bracketed difference between the geometric factors is always non-negative (e.g., adding more sensors cannot reduce the source screening probability), so that: igure 7: Optimal and sub-optimal azimuthal solutions for extra sensor locations to supplement an existing sensor deployment. (a): Solution curves to eq. 28 that show the dependence of optimization parameter Λ on the geometric deployment factor (left factor in eq. 21; vertical axis) on azimuth ψ (horizontal axis). Points below the dashed, "stationary" line marks azimuthal regions where both ∆Λ and ∆η increase. The global minima (circles) mark solutions that maximize source screening parameter Λ. The local minima (triangles) mark solutions for the best location of an extra receiver that is near the worst deployment location (squares). (b): Screening probability curves Pr D (vertical axis) for the decision rule in eq. 28 when deployments include one and two extra sensors to supplement an existing N = 16 station network (gray markers, upper left), compared against the scaled faulting signal (horizontal axis). The radiation pattern source has the same focal mechanism as the 2003 Lop Nor tectonic event (beach ball at top right). The thick gray curve labeled with a gray, circled 0 marks the performance of the screening statistic that associates with the original instrument deployment (gray circles). The other numbered curves show the screening performance of the optimally deployed N = 17 and N = 18 networks. These curves associate with the sensor locations that we depict with deployment schemes (upper left) and minima in (a). Squares mark curves that show worse-case deployment azimuths that correspond to maxima on the curve in (a).
The last term of eq. 28 defines the objective function for Rayleigh wave sampling. It omits explicit dependence on data R. Therefore, the optimal strike-relative deployment azimuth estimateψ +1 for a receiver when σ 2 M is unknown depends only on current receiver placement and the relative seismic wave structure (α and β). Other values of sensor placement where A H T +1 H +1 −1 A T < 1 increase Λ from its prior value (eq. 21). Auxiliary sensor azimuths where A H T +1 H +1 −1 A T > 1 decrease Λ from its prior value. To quantify the change in screening power, we reapply eq. 23 to estimate the updated threshold η and reapply eq. 22 to compute source screening rates with N + 1 sensors; we then update the decision rule (eq. 18). Fig. 7a shows π-periodic optimal and sub-optimal solutions for extra receiver locations that supplement 16 existing receivers uniformly deployed along a 75 • azimuthal arc (gray circles, Fig. 7b). The objective function shows local and global extrema of The maxima occur at deployment azimuths that minimize Λ and therefore the probability of screening binary source types; these sensor locations provide the poorest improvement in source-type screening. Local minima (triangles) provide a locally optimal sensor placements near these worst-case sensor locations. These local minima solutions may apply to scenarios in which deployment outside restricted azimuthal ranges is impossible (e.g., inaccessible areas). Our global solutions (circles) confirm expectation that optimal receiver deployment should include points that include no current receiver coverage. These optimal deployment locations are not orthogonal to the mid-point of the azimuthal arc, but depend on the relative values of the body wave speeds (parameter c). Fig. 7b shows screening curves for the decision rule 18 when a receiver network contains the backbone of existing sensors (inset gray circles), supplemented with poorly placed sensors (squared one and two), versus optimally placed extra sensors (circled 1 and 2). The optimally placed sensor that supplements this 16 sensor backbone network forms a 17 sensor network the improves screening rates of the faulting source by > 6× that of the backbone network, for sources with the greatest faulting signal shown (circled 1). An additional, optimally placed sensor (circled 2) provides a greater gain in screening performance, relative to an 18 sensor network that includes two poorly placed auxiliary sensors (squared 2). The screening performance our test achieves with this azimuthally limited, 18 sensor arrangement outperforms average test that uses the 24 sensor deployment in Fig. 6.
Lastly, we emphasize that the update to η maintains a consistent false-attribution probability for the discriminant before and after adding any additional receivers. Failure to update this threshold when ∆Λ < 1 could lead an observer to erroneously conclude that adding receivers at particular values of ψ +1 (shaded domain of Fig. 7(a)) can reduce the screening power of the same test.

Case I application: threshold faulting moments and missed detections
We now quantify the largest faulting source signal (as defined by (DS + SS) 2 /σ 2 M ) that could be mistakenly attributed to an explosion, at a fixed false attribution probability Pr F A . This application is particular important to underground explosion monitoring applications because it provides fundamental physical limits on the the discrimination power of eq. 16 tests to screen explosive from non-VDS faulting sources.
To compute this signal, we first note that the probability of mistaking a non-VDS faulting source for an explosive source is equivalent to the event that L θ (R) < η when H 1 is true (eq. 18). Here, threshold η relates to the CDF under hypothesis H 0 through eq. 23. The probability of miss-attributing a non-VDS faulting source to an explosion then equates to 1 − Pr D and relates to the faulting signal through a nonlinear equation that requires inversion of the noncentral-F CDF for Λ: The term 1 − Pr D is equivalent to a "miss" in waveform detection theory. The faulting signal SNR that solves eq. 29 relates to noncentrality parameter estimateΛ through: Our estimate for the source term ofΛ, as measured by a given sensor deployment locations ψ m (m = 1, 2, · · · , N ) that are stored in H, is:  Moreover, the variability in such threshold estimates is asymmetric: the log-scale maximum threshold estimates (top curve) depart from the mean (middle curve) significantly more than do log-scale minimum threshold estimates (bottom curve). Noisy radiation patterns that associate with these solutions in Fig. 9(b) show that the decision defined by eq. 18 can mistake (with probability 1 − Pr D = 0.1) the subtle four-lobed shape for Figure 8: Three solutions curves to eq. 31 that quantify the largest strike-slip faulting source signal that the decision rule in eq. 18 can mis-identify to originate from an explosion with probability 1 − Pr D = 0.1. In each case, N = 12 sensors confined to three distinct azimuthal deployment ranges (at right) sample the strike-slip Rayleigh wave radiation pattern (beach ball inset). Each curve compares faulting signal SNR (vertical axis) against false attribution probability (horizontal axis), and error bars measure standard error from 100 random sensor deployments. Filled circles (right) mark a particular, random deployment (we use standard error here because standard deviation markers reduce readability). Three noisy radiation patterns of a strike-slip source (eq. 16) with three corresponding azimuthal gaps in sensor sampling. The decision rule in eq. 18 samples these radiation patterns at N = 12 locations to falsely attribute the source of each pattern to an explosion, with probability 1 − Pr D = 0.1. A four-lobe Rayleigh wave radiation pattern is evident in the case that sensors sample the radiation field without azimuthal constraint (blue pattern). We emphasize that the decision rule in eq. 18 mistakes this pattern as circular (with a 0.1 probability).
a noisy circular radiation pattern, even with full azimuthal sensor coverage.

Case II: σ R and θ unknown; σ M known
We progress from Case I and now assume that a sensor network records a Rayleigh wave radiation pattern with unknown random process variability (σ 2 R in eq. 16), unknown faulting parameters (θ in eq. 11), and a stochastic variance σ 2 M that is either known, or has a significantly lower estimator variance than than that ofσ 2 R . This case applies to idealized monitoring scenarios in which an underground testing complex colocates with a shallow fault of unknown geometry and radiation variability, but with well-characterized stochastic variability σ 2 M that is perhaps quantified from historical records of proximal explosions. Sensor networks that record a low magnitude seismic event that originates from this test site faces the same challenge as that in Case I, but the screening statistic has a distinct form. We treat Case II in less detail than Case I because it appears less practically applicable than Case I while remaining more complicated. Appendix C derives several results that we state here. For Case II, the log-GLRT equivalent to eq.
where eq. 14 and eq. 15 define f R (R; H 0 ) and eq. 13 defines f R (R; H 1 ). The MLE in which the text that follows eq. 19 defines the projector matrix P ⊥ H . Our eq. B.2 defines the MLEsθ k for parameter vectors θ k (k = 0, 1). We input these parameters to eq. 33 and reuse the symbol L θ (R) to write the resultant, equivalent GLRT statistic (algebra omitted): We use the MLEs forθ k to express L θ (R) as a sum of two statistically independent random variables (X and Y ): We interpret the first quadratic form ( P ⊥ H R 2 //σ 2 M ) in eq. 35 as the residual radiation pattern energy not attributable to the Rayleigh wave model, relative to random noise variance. The second term ( P X R 2 /σ 2 M ) is energy due to non-VDS faulting, relative to the same noise variance. The decision rule that tests the function of these two terms Figure 10: Screening curves for Case II (thick, solid curves) compared with Case I test statistics (thin, dashed curves). The plot format is similar to Fig. 5: Case I and Case II curves present 100%× the probability Pr D (vertical axis) that eq. 18 and will screen earthquakes with three distinct focal mechanisms from an explosion, using N = 12 sensors, all plotted against faulting signal SNIR (horizontal axis; eq. 21). Each curve is an average of 100 screening curves (eq. 24) that associate with sensor deployments that uniformly sample azimuths with a 90 • gap along the dashed circle centered at the hypocenter of the shallow source (orange star); the filled circles mark a particular, random deployment. The topmost curve (blue) presents such an average Case II screening curve for an earthquake with a focal mechanism that matches that of "Event 1" from the Rock Valley earthquake; the second curve (red) presents a similar average for an earthquake with a focal mechanism like that of the Kara Sea event; and the bottom curve presents a similar average for an earthquake with a focal mechanism like that of the 2003 Lop Nor event. The shaded region indicates ± one standard deviation (±σ) from the top curve mean. All curves correspond to Pr F A = 10 −3 . in eq. 35 is: The thresholdη in eq. 36 is strictly an estimate that maintains a false-attribution probability Pr F A , consistent with Case I. We describe our estimation strategy for this threshold in Appendix C (see eq. C.5). To compute the screening performance for Case II (SNIR versus Pr D ), we note that the terms X and Y are independent, piecewise invertible functions. In particular, the term X is piecewise invertible in argument s = P ⊥ H R 2 /σ 2 M since function x(s) = s − N ln(s) monotonically decreases over s < N and monotonically increases when s > N . The ratio s ∼ χ 2 N −3 (0), where χ 2 N −3 (0) is chi-squared distribution function with N − 3 degrees of freedom and a zero noncentrality parameter; the distributional form of s is identical for both circular and non-circular radiation patterns. Unlike the Case I test statistic, these PDFs are well defined if N ≥ 4 sensors rather than eight sensors. The ratio Y ∼ χ 2 1 (Λ) is a chi-squared distribution function with 1 degrees of freedom and a nonzero noncentrality (Λ > 0) under H 1 .
We use these results to compute the PDF of the sum Z = X + Y through a variable transformation on the original, statistically independent random variables. Given that X and Y respectively have PDFs f X (x) and f Y (y), the PDF f Z (z; H i ) for Z is their convolution f X (x) * f Y (y) under hypothesis H i (i = 0, 1). Omitting the explicit algebra, we symbolize the PDF f Z (z; H i ) as: We compute eq. 37 in the frequency domain with the characteristic function method (eq. C.4). The PDF f χ 2 1 (•) includes a finite noncentrality parameter under the alternative hypothesis (i = 1). The theory we detail in Appendix C demonstrates that this parameter equates to the same scalar Λ in eq. 21 that we derived from Case I. This means that Λ completely quantifies the screening capability of L θ (R) to test between H 0 and H 1 .
The PDF f X (x) is identical under H 0 and H 1 . We interpret its effect on f Z (z; H i ) as a smoothing window in the convolution operation in eq. 37 that is independent of the Rayleigh wave radiation pattern model. Fig. 10 shows the screening power of the Case II test statistic compared against that of the Case I test statistic, for three source types and a fixed deployment constraint (an azimuthal gap of 90 • ) for N = 12 sensors. The graph format is identical to that shown by Fig. 4, Fig. 5 and Fig. 6, in which each curve represents an average (eq. 24) over 100 deployment realizations. Each such realization effectively recomputes the deployment factor A H T H −1 A T −1 of Λ and re-convolves the PDFs f Y (•) and f X (•) in the frequency domain as a spectral product. We document the additional numerical details that we implemented to construct the screening curves shown by Fig. 10 in text following eq. C.5.
Importantly, these curves demonstrate that the Case II decision rule provides a superior source-screening capability, at least for the same faulting signal SNIR and thresholds required to maintain a false attribution probability of Pr F A = 10 −3 . We expect the Case II test statistic to provide such an improved, average performance because the Case I statistic additionally requires an MLE of the noise variance (σ 2 M ) under H 0 ; in general, test statistics with unknown parameters show poorer screening performance than competing test statistics with fewer unknown parameters (Kay, 1998, pg. 195). The variability in screening performance with sensor placement over permissible azimuths (note the 90 • gap), however, is comparable between both Case I and Case II (note shading). Both tests are also similar in the relative order of screening rates with focal mechanism. In particular, the Case II statistic shows a qualitatively similar ordering in screening probability when compared to the Case I statistic. The Case II decision rule (eq. 36) that tests radiation pattern data that resembles the Rock Valley earthquake (for example) shows a much higher probability of correctly screening such events from isotropic sources than, say, data that records an event that resembles the 2003-03-13 Lop Nor earthquake.

Discussion
Our binary hypothesis test (eq. 16) quantifies an observer's idealized ability to screen shallowly buried explosion sources from non-VDS faulting sources. These tests effectively measure deviations from circularity in Rayleigh wave radiation patterns that any faulting terms induce at the source location (Fig. 1). While the binary test provides only a commensurate, binary decision on the significance of any non-circular contributions to a measured radiation pattern, it affords significant insight into an observer's capability to distinguish source types.

Case I
We first consider our results under the general assumptions of Case I, when σ 2 R = 0: • An observer can use azimuthally distributed receivers that sample a noisy Rayleigh wave radiation pattern to form a test statistic and decision rule (eq. 18). This rule screens explosive from non-VDS faulting sources at a given threshold η that maintains a fixed false attribution probability Pr F A .
• We interpret this screening statistic as a ratio between radiation pattern energy due to non-VDS faulting and residual radiation pattern energy that is not captured by our Rayleigh wave model. This test statistic has a noncentral F -distribution that quantifies this observer's screening capability, and is well-defined only when an observer deploys eight or more receivers.
• Our test between circular versus non-circular Rayleigh radiation patterns is equivalent to a general Gaussian signal detection problem (eq. 18). The product of a faulting SNR term and a scalar deployment term (eq. 21) defines the test screening power. The first term measures the squared sum of the dip-and strike-slip faulting components, relative to moment noise: (DS + SS) 2 /σ 2 M . The deployment term A H T H −1 A T −1 depends on azimuthal sensor distribution and the body wave speeds of the host medium (parameter c).
• Strike-slip sources provide the highest faulting signal strength M 2 0 /σ 2 M for a given faulting moment M 0 . Eq. 18 therefore screens strike-slip from explosive and vertical dip-slip faulting sources with the highest probability Pr D = 1−F 1,N −3 η, M 2 0 /σ 2 M , among other sources with the same faulting moment M 0 and constant false attri- • Radiation pattern screening probabilities depend on three parameters: azimuthal gap, source radiation pattern, and the number of deployed sensors. Fig. 4, Fig. 5 and Fig. 6 show that certain combinations of these three parameter sets produce qualitatively similar radiation pattern screening curves.
• An observer can use the screening test to estimate the azimuthal placement of additional sensors that supplement an existing network, and which maximize the probability of screening circular from non-circular radiation patterns. Moreover, deploying additional sensors at any azimuth can never reduce screening rates.
• Lastly, observers have a non-negligible probability of mis-attributing non-circular radiation patterns of faulting sources to explosions when they measure these patterns with sparse sensor deployments. The SNR of such a faulting signal increases exponentially as the azimuthally gap in sensor deployment exceeds ∼ 90 • . This SNR can exceed an order of magnitude as gaps increase from 0 • to 135 • (Fig. 9).
We next consider our (more limited) results under the general assumptions of Case II.

Case II
The assumptions under Case II explicitly hypothesize that an unknown, random process superimposes with a deterministic radiation pattern signal and imprints on the total pattern. It also assumes that stochastic variability within the radiation pattern (noise) is effectively known under each hypothesis. This latter assumption remains less practical because the seismic noise field is often insufficiently stationary to be well-characterized, at least at multiple azimuths. We can, however, report several outcomes: • The Case II statistic increases with Rayleigh wave energy sourced by non-VDS faulting, relative to noise power ( P X R 2 /σ 2 M ). It is non-monotonic in residual energy that is not captured by our Rayleigh wave model and normalized by noise • The screening power of the Case II decision rule (eq. 36) depends on a scalar function of faulting signal and sensor deployment. This scalar term Λ is identical to that under Case I (eq. 21).
• A spectral product of two PDFs quantifies the screening performance of the Case II test statistic. The resulting PDF does not appear to have a known standard form. It's CDF predicts that the Case II decision rule (eq. 36) is more effective at screening circular from non-circular radiation patterns than the Case I decision rule (eq. 18), because it requires fewer MLEs.
• The Case II test statistic is (as opposed to Case I) well defined with N ≥ 4 sensors (rather than eight sensors).
Both the Case I and Case II results have implications to screen SPE-triggered explosionplus-damage signals from historical faulting signals that are each sourced in the NNSS testing complex. We suggest that an observer can use N = 12 sensors deployed around the source with an azimuthal gap of 90 • or less to screen faulting signals that match such historically recorded mechanisms, provided these data have sufficient SNR (Fig. 10). A confirmation of screening power from such a validation exercise remains necessary if this screening statistic holds efficacy in non-idealized settings.

Caveats
Our results remain subject to at least four more specific caveats. First, our work applied a planar model applicable to local distances that concern records from sources like the SPE or New Hampshire shot series. More established observational and theoretical work that quantifies the imprint of faulting motion on Rayleigh wave radiation patterns (e.g., tectonic release) considers far-regional to teleseismic data. Such Rayleigh waveforms are dominated by long period spectral content and propagation distances of many hundreds to thousands of km. We assert hypothesis tests with our model remain applicable to such global problems. Specifically, we argue that we can rotate the data at each sensor location through multiplication by a unitary matrix to a locally planar coordinate system; this idempotent matrix multiplication does not change the distributional form of such quadratic forms, nor the expected performance of each test statistic.
Second, we did not explore the relative contributions of the explosion versus the CLVD source to the circular part of the radiation pattern in our model. This CLVD term can exceed the explosion term in size, which differs in sign so that Rayleigh waveforms reverse polarity. While the Case I and Case II test statistics exclude any such circular terms (Λ excludes ∆R), we concede that a zero-mean damage signal may appear in the data as a random process. Physically, this means dilatation that is symmetric about the vertical axis may manifest as a late time CLVD-source. Our linear model does not accommodate such moment-dependent changes in the total data variance. The possible presence of such a CLVD signal motivates our future consideration of the Case II test statistic.
Third, we assumed that noise recorded at distinct sensor azimuths is uncorrelated, largely for mathematical convenience. Correlated noise reduces the effective degree of freedom that parameterize the distribution of each test statistic. This means that our theoretical screening rates will exceed our observed detection rate if our method does not accommodate for noise correlation. We can make such accommodations by computing empirically-validated estimates for the effective degrees of freedom in observed data (Carr et al., 2020, eq. A1, eq. A3). This argument applies to both Case I and Case II.
Fourth and lastly, our noisy Rayleigh wave models impart "ragged edges" to the graphical radiation pattern representations (Fig. 9). These distortions from their deterministic patterns that are sourced by noise may be relatively small when compared to distortions sourced by deterministic focusing or other scattering effects not included in this model (in a real Earth). Such deterministic distortion would shape the radiation pattern from isotropic source (a pure explosion) so that this pattern would appear to have a partial "lobe". Our decision rules (eq. 18 or eq. 36) would then show increased Type 1 errors, that is, a higher probability of falsely deciding that the data include a non-VDS faulting signal (choosing H 1 when H 0 is true). Alternatively, scattering of Rayleigh wave energy from surface topography that is comparable to Rayleigh wavelengths might homogenize the radiation pattern from a source with a considerable faulting signal. We expect this homogenization far from the source through energy conservation and divergence arguments: more energy is available to be directed from anti-nodal azimuths to nodal azimuths from any non-zero scattering effects. This example would elevate Type 2 error rates, that is, the probability that our decision rules incorrectly decide that a tectonic source is explosive in origin (choosing H 0 when H 1 is true); this latter example would imply that eq. 18 will mistake sources with larger faulting signals than shown in Fig.   9 as explosions more often than we predict. In the language of signal detection, these collective physical processes present in a real Earth will inflate false alarm and missed detection rates over our predicted rates.

A Physical Source of a Screening Error
We explicitly treat an analytical example of an erroneous screening test that we apply to a radiation pattern that is induced by a vertical crack (see comments in Section 4.2) opening within a stratified half-space. This example produces a noncircular radiation pattern without the faulting signal that is modeled by hypothesis H 1 . To construct this pattern, we first consider the moment density tensor for a shallow displacement discontinuity (Pujol, 2003, eq. 10.6.6.). We then write the (full) moment tensor for such a discontinuity that is source by a vertically oriented surface crack of area A, located at the source point ξ 0 , that opens in tension a distance u(ξ 0 ) in the Cartesian y direction (algebra omitted): Eq. 38 omits the frequency-domain representation for the source time function that would describes the history of the crack-face displacement and appear as a factor of M . Figure 11: A range of predicted screening curves (p-box) superimposed with the mean observed screening curve for the Case I decision rule (eq. 18), which we apply to a synthetic radiation pattern collected from azimuthal records of Rayleigh waves triggered by an opening surface crack (beachball at center, right; see Section 5.4). Curves present 100%× the probability Pr D (vertical axis) that eq. 18 will use N = 24 sensors to decide that a fracture source produces a noncircular radiation pattern and includes a faulting signal, plotted against radiation pattern SNR (horizontal axis). The shaded region shows a measure of the range of theoretical screening curves. The curve that marks the lower limit of this range is parameterized by noncentrality parameter Λ = P X R d 2 /σ 2 M − σ, where σ is the standard deviation marking expected variability of estimates for Λ, which originate from σ 2 M . The curve that marks the upper limit of this range is similarly parameterized by noncentrality parameter Λ = P X R d 2 /σ 2 M + σ. The dark curve presents an average of 2 × 10 4 synthetic decision rule counts, that is, the number of times that eq 18 chooses hypothesis H 1 . All curves correspond to Pr F A = 10 −3 . The synthetic "ragged edged" radiation patterns with variance σ 2 M at bottom associate with the labeled SNR values. White circles mark azimuthal sensor locations.
We note that this moment tensor is a sum of a tensor from an explosion source, and a tensor from a force couple that is aligned with the opening crack. We use the M in eq. 38 and the radiation pattern that appears in both eq. 2 and eq. 3 to write the Rayleigh surface displacement u as (again, omitting the source-time function): The angle ϕ in Eq. 39 measures azimuth from the linear crack face. We pair each term with the deterministic radiation pattern coefficients that associate with trigonometric basis functions in the noisy description of R in eq. 9: which have units of moment (F · m). The term ρA u(ξ 0 ) β 2 represents the elastic energy released by volumetric crack growth. Such sources are common in seismogenic hydrofracture events within glaciers (Carr et al., 2020;Taylor-Offord et al., 2019) and ice sheets (Carmichael et al., 2015).
We now suppose N sensors distributed at distinct azimuths along a circle that is centered at the source crack samples the radiation pattern in eq. 40, which we store in vector R d (subscript d indicates deterministic). These sensor azimuths populate the system matrix H in eq. 19, and the noncentrality parameter equivalent to eq. B.5 becomes P X R d 2 /σ 2 M ; here, eq. 19 defines P X and σ 2 M is Gaussian noise variance. We then add noise with this variance to the fracture-generated radiation pattern and process these synthetic data (now R = R d + n, where n is noise) with our Case I test statistic L θ (R) and decision rule (eq. 18) over a grid of crack release energies and SNR values.
In particular, we considered ratios of 20 ≤ ||R d || 2 /σ 2 M ≤ 400 and N = 24 sensors. Fig.  11 shows that the Case I test erroneously screens opening cracks as explosions for low SNR, or as non-VDS faulting sources at higher SNR values. Either decision may be interpreted as an error.

Conclusions and Future Work
The goal of this work was to determine if a discriminant for shallow source types that tests Rayleigh wave radiation pattern shapes is justified on theoretical grounds. We conclude that this approach is likely not practical in most passive monitoring scenarios that require sampling the Rayleigh wave field over a non-horizontally stratified medium with heterogeneities and substantial topographic relief (a real Earth) for two reasons.
First, an observer maintains a good discrimination capability between circular and noncircular radiation patterns only with a large number of sensors with good azimuthal coverage and high SNR in a cylindrically symmetric medium. Second, any distortion to the Rayleigh radiation pattern sourced by deterministic processes (not just noise) will increase Type 1and Type 2 error rates in the decision rules that provide a source-type discrimination capability, over that of the ideal cases.
In particular, our hypothesis tests on radiation pattern shape show our tests achieve good, but idealized screening power over only certain regions of parameter space that include deployment azimuth, focal mechanism, and sensor number. We demonstrated that a Rayleigh wave radiation pattern test that exploits a sufficient number of sensors (> 12) with limited azimuthal gap (< 90 • ) to record radiation from sources with moderate strike-slip faulting signals (SNR > 20), like the historical Rock Valley Event 1, show a high probability Pr D of success (Pr D > 0.9). The probability of such success diminishes with more ambiguous focal mechanisms that resemble certain historical, shallow earthquakes that locate near the Lop Nor nuclear test site in China. Our further technical developments require validating the screening capability of this method against SPE shots that source near historical, tectonic sources of non-circular Rayleigh wave radiation.
To overcome challenges of limited explosion data (single shots), we will perform semiempirical tests on existing explosion records. These tests will infuse amplitude-scaled Rayleigh waveforms triggered by a known source into long noise records, thousands of times, and will re-compute screening curves from the resultant, semi-synthetic data. This process thereby will provide a set of "observed" screening curves that we can directly compare to our predicted screening curves that test radiation pattern circularity.
In achieving the goal of this work, our hypothesis testing approach did provide two unexpected, ancillary benefits to experimental planning of seismic deployments. First, we can quantify how a particular geographical placement of additional sensors to supplement an existing dense network can drastically increase our capability to screen the non-VDS faults from the explosive sources. We conclude from this case that optimal sensor deployment does not depend on data records, only prior sensor locations and the relative compressional versus shear wave speed. Second, our screening statistic quantifies the largest faulting signal that an observer could expect to record and misattribute to an explosion source (for some probability). The size of this source increases rapidly with azimuthal gap in receiver deployments that exceed 90 • . Our estimates are applicable to test-ban treaty verification scenarios that require attributing discriminant evidence to source identity, but (as stated) we must first assess if such a test even has marginal screening power with real data.
We concede that there are at least three technical considerations outside the scope of this current study (additional to caveats). First, we did not supplement our radiation pattern tests with parallel analysis of Love wave radiation pattern circularity (present under hypothesis H 1 ). Second, we do not consider far-regional and teleseismic observations that require spherical propagation models. Third, we did not explore how damage processes that contribute to the CLVD source, which can reverse Rayleigh waveform polarity, contribute to the apparent stochastic noise in our data. We note that some evidence suggests that random processes can model vertical axis, cylindrically symmetric damage signals in Rayleigh waveform records. Such processes will contribute to our Case II radiation pattern models.
While this theoretical analysis did not use not waveform data, our results do suggest that radiation pattern shapes provide screening power as a physically-interpretable, sourcetype discriminant in some idealized cases. Such cases include controlled, laboratory experiments with engineered structures, or natural, horizontally stratified structures that are leveraged in seismic studies elsewhere (Biagi et al., 1990;Hudson et al., 2020;Knox et al., 2016;Lott et al., 2020;Lu et al., 2007;Toksöz et al., 1971). If future application of our tests with locally recorded data prove successful, we will begin testing it on regional, low frequency data. If our method instead fails, then we will quantify the geophysical source of this disagreement to gain physical insight into how Rayleigh radiation pattern shape can or cannot be tested as a discriminant.

Data and Resources
Our theoretical study used no data. shows Rayleigh wave radiation patterns (black curves) with their RMS values R (purple curves), superimposed withR ± σ R (red, dashed lines). The total seismic moment (M 0 + M I ) scales all radiation pattern plots. We emphasize that figure labels σ R in this figure exclude (due to lack of visibility), but differs from the variance present in an unknown random process, as considered under our Case II analyses.

A Deterministic Radiation Pattern Moments
We compute deterministic moments of a Rayleigh wave radiation pattern R that is produced from buried, shallow sources; we emphasize that "moments" in this appendix often refer to probabilistic moments, so we explicitly label seismic moments. The probabilistic moments include the (1) deterministic meanR and (2) the deterministic variance σ 2 R of R. We also bound their behavior with their root-mean-square (RMS) values R 2 , and σ 2 R to measure the expected deviation of the radiation pattern fromR + σ 2 R . We assume sensors are deployed with equal probability around the source so that ψ is uniformly distributed over [0, 2π].
To compute both moments, we assume the source is embedded within a vertically stratified half-space and surrounded by an imaginary, unit radii cylinder to a depth equal to the dominant wavelength λ of the radiated elastic energy; the current context provides little risk to confusing wavelength with screening parameter Λ and we therefore reuse this conventional symbol. Regardless of symbolism, we then integrate the squared radiation pattern R 2 over our cylinder, and normalize the integrand over the cylinder surface area. This process is analogous to using Gaussian surfaces in electrostatics to compute electric fields (Griffiths, 2017, Chapter 2). In our case, this integrated radiation pattern is (from eq. 7): The radiation coefficient weights [1, cos(2ψ), sin(2ψ)] in eq. A.1 are mutually orthogonal over the azimuthal interval [0, 2π]. Hence, cross terms such as DS cos(2ψ) · SS sin (2ψ) integrate to zero. The remaining terms reduce eq. A.1 to: Comparing to eq. A.1, we conclude: The covariance σ R i R j = cov (R i , R j ) similarly quantifies how samples of the deterministic radiation pattern correlate at distinct azimuths ψ i and ψ j : Products like cos(2ψ i ) sin(2ψ i ) also integrate to zero. Therefore, for i = j: Radiation pattern observations collected from distinct azimuths are therefore statistically independent, when the prior distributions on ψ i and ψ j are uniform. Because this model is a mathematical idealization, we concede that observations collected from sufficiently proximal azimuths are likely to be correlated. We can conclude, however, that the total covariance in an observed radiation field is representable as σ 2 M I when sensor deployments are sufficiently sparse. Fig. A.1 illustrates relationships between deterministic radiation pattern moments for several fault types.

B Derivation of the Case I Screening Statistic
We apply the hypothesis test of eq. 16 and combine eq. 13 through eq. 17 to compute an equivalent GLRT L The MLEs for the unknown parameters are (Kay, 1993, pg. 252, eq. 8.52): in which i = 0, 1 indicates distinct noise variance estimates under hypothesis i, and σ 2 R = 0 when i = 0. To reduce L (2) θ (R) into a more interpretable test statistic, we rewrite the numerator and denominator of eq. B.1 using projector matrices. These matrices are documented in eq. 19 and the text that immediately follows eq. 19. With these definitions, the Case I test statistic is a ratio of two statistically independent quadratic forms: so that P X projects onto rank-1 subspace span{X}. Matrices P ⊥ H and P X project onto orthogonal spaces since P H X = X. The Pythagorean identity therefore implies that P ⊥ H R + P X R 2 = P ⊥ H R 2 + P X R 2 . This factorization simplifies the screening statistic into a reduced form L θ (R) that is ratio of two statistically independent terms.
We then scale L θ (R) by the ratio of the alternative hypothesis to null hypothesis signal variance (σ 2 M + σ 2 R )/σ 2 M : The test statistic L θ (R) has noncentral-F distribution parameterized by scalar Λ when R k includes additive Gaussian noise (1 ≤ k ≤ N ). We write its distributional dependence as L θ (R) ∼ F 1,N −3 (Λ > 0) and the cumulative, noncentral F distribution for L θ (R) as eq. 20. Scalar Λ has the general quadratic form when Aθ = b under H 0 : The scalar Λ completely quantifies the screening capability of L θ (R) to test between H 0 and H 1 . Specifically, increasing values of Λ decrease the overlap between the density functions for test statistics that describe cylindrically symmetric (Λ = 0) and faulting sources (Λ > 0). To factor Λ into its source-dependent and deployment depend parts, we distribute the scalar A H T H −1 A T from (DS + SS) 2 . This factorization results in eq. 21.

C The Density function for Case II
This appendix details the development and performance of the Case II test statistic L θ (R) (we redundantly label the Case I statistic with the same symbol). We first note that the multi-valued inverse of x(s) that we write as s(x) in the PDF arguments of eq. 37 is not a proper function. Rather, x(s) = s − N ln(s) is minimized at s = N , monotonically decreases for all s < N , and increases for s > N . The Lambert function W (•) compactly represents this multi-valued inverse s(x): in which the properties of W (•) constrain s(x) to be non-negative. To implement the change-of-variables in eq. 37, we write the PDF f X (x; H i ) as a two-term sum over each single valued branch of the Lambert function (i = 0, 1): estimates for PDFs under the competing hypotheses. We selected Λ = 10 under H 1 and a false attribution probability of Pr F A = 5 · 10 −3 . The FFT inversion appears stable when Λ is not too large (Carmichael et al., 2020, Section S2). Sources that trigger radiation patterns with high SNR faulting signals will be obvious to screen, and we therefore do not consider this restriction on Λ of practical significance. To compute screening probability values for each of the focal mechanisms shown in Fig. 10, we integrate the PDFs in eq.
C.5 over a 512 point grid to construct a CDF for each source focal mechanism. We then repeated this integration over a 500 point grid of faulting SNIR values (horizontal axis).
We then estimate each screening probability by interpolating the CDF at a screening threshold that is consistent with a false alarm value of Pr F A = 10 −3 .