Likely PKS-PKP from Array Processing of Noise Records

Seismic noise has been widely used to image Earth’s structure in the past decades as a powerful supplement to earthquake signals. Although the seismic noise ﬁeld contains both surface-wave and body-wave components, most previous studies have focused on surface waves due to their large amplitudes. Here, we use array analyses to identify body-wave noise traveling as PKP waves. We ﬁnd that by cross-correlating the array-stacked horizontal-and vertical-component data in the time windows containing the PKP noise signals, we extract a phase likely representing PKS-PKP, the diﬀerential phase between PKS and PKP. This phase can potentially be used for shear-wave-splitting analysis. Our results also suggest that the sources of body-wave noise are extremely heterogeneous in both space and time, which should be accounted for in future studies


Introduction
Recent decades saw a rapid expansion of studies using seismic noise to image Earth structure (e.g.Shapiro et al. (2005), Bensen et al. (2007), Brenguier et al. (2008), Lin et al. (2009), Poli et al. (2012), Nakata et al. (2015)).Most of these studies focused on extracting surface-wave signals from the noise field because surface waves usually dom-inate the signals retrieved by noise cross-correlation.This observation is commonly attributed to the prominence of surface waves in Earth's noise field as a result of noise sources, such as wind and ocean waves, occurring mostly at the surface.Despite their lower amplitudes, body-wave signals have occasionally been retrieved from noise cross-correlations and used to image Earth structure (e.g.Poli et al. (2012), Nakata et al. (2015), Feng et al. (2021)).A major advantage of body waves over surface waves in studying Earth structure is that body-wave reflected and converted phases are sensitive to material discontinuities in Earth's interior (e.g., the Moho and the core-mantle boundary (CMB)), which cannot be resolved with surface-wave data alone.However, body-wave reflection and conversion signals are weaker than direct phases and thus more difficult to observe in the cross-correlation functions, which are typically nosier than earthquake recordings.Therefore, techniques capable of enhancing body-wave reflection and conversion signals are needed to better image Earth's discontinuities with noise records.
In addition to imaging using seismic noise, in recent years major advances have been made in understanding the sources of Earth's noise field (e.g., Gualtieri et al. (2014), Nishida and Takagi (2016), Liu et al. (2020), Retailleau and Gualtieri (2021)).Many contributions were made by studying body-wave noise signals with array techniques (e.g., beamforming and back-projection), which suggests that weak body-wave noise signals can be enhanced with array processing to better image Earth structure.These studies also showed that body-wave noise sources, which are usually associated with storms in the oceans, are likely spatially and temporally heterogeneous, which implies that body-wave signals could be better retrieved through seismic interferometry if the variations of the bodywave noise sources are properly accounted for.
Here, we present observations of body-wave noise propagating as PKP using data collected by a dense broadband seismic array in the central US.We further show that a phase likely representing PKS-PKP can be extracted by cross-correlating the arraystacked horizontal-and vertical-component noise records in the time windows containing the PKP noise signals.We then discuss the potential applications of this seismic phase and the implications of our findings for seismic interferometry.

Data and preprocessing
We mainly use the continuous data collected by the Ozark Illinois Indiana Kentucky (OIINK) Flexible Array Experiment (network code: XO), a dense 2D broadband seismic array with a station spacing of ∼ 25km located in the central US (Fig. 1).To make the resolution of our results more isotropic, we select the OIINK stations located in a 100-km radius circle and also include the Transportable-Array stations in this range (Fig. 1b).Although the two arrays together span 2011-2015, to ensure a reasonable resolution we focused on time windows with more than 20 active stations, which limits our analysis to a roughly one-year period between June 2012 and August 2013.We downloaded the continuous data from the IRIS Data Management Center in one-hour time time windows, removed the instrument response, and band-pass filtered the data to 2-10 seconds, which contains the secondary microseism energy (Retailleau & Gualtieri, 2021).
To avoid the effects of earthquakes and instrument malfunctions, we removed the 1-hour time windows containing the first arrivals of global earthquakes with magnitude > 5 and those containing amplitudes >1 × 10 −5 m s −1 .

PKP signals from beamforming analysis
We performed conventional linear beamforming with all three components (vertical, east, and north) of our array data to characterize the directional properties of the noise field.To save computational cost, we first performed a reconnaissance analysis over the slowness range ±0.2 s km −1 at a grid spacing of 0.013 s km −1 in the W-E and S-N directions.The resulting vertical-component slowness images clearly show beams with slowness <0.04 s km −1 , which likely represent PKP signals (Fig. 2a).The horizontal-component slowness images also show local maxima corresponding to the PKP beams on the vertical component, though the background noise is significantly higher on the horizontalcomponent images (Fig. 2a), which could be due to the near-vertical particle motion of PKP or a more homogeneous distribution of horizontal-component noise sources.The slowness images of some time windows also show multiple peaks (e.g., 2013-07-06-00-00-00; Fig. 2a).
For seismic imaging, we prefer to use time windows dominated by PKP energy from a single direction because this resembles that of earthquake sources, which may make techniques in earthquake imaging readily applicable.To identify these time windows, we find the maximum in the PKP range (slowness <0.04 s km −1 ) of each vertical-component slowness image and the corresponding slowness vector, which we refer to as the PKP slowness.We then define the vertical-component normalized PKP -beam amplitude as the ratio between the maximum amplitude in the PKP range and the average amplitude of the whole slowness image, which measures the power of the strongest PKP beam relative to the background noise.We further define the corresponding normalized PKPbeam amplitudes for the horizontal components as the ratios between the amplitudes at the PKP slowness and the average amplitudes of the whole slowness images.We finally define the three-component normalized PKP -beam amplitude (hereafter "PKPbeam amplitude") as the product of the normalized PKP -beam amplitudes for the three components.We regard the time windows with PKP -beam amplitude > 2, which account for about 10% of all the time windows, as windows dominated by PKP energy from a single direction and make a histogram of the PKP slowness of these time windows, which shows that the vast majority of these time windows have slownesses close to the b and c caustics of PKP (Fig. 2b).This phenomenon is probably due to the amplification of PKP near the caustics.
To identify the source locations of these PKP beams, we performed beamforming for the vertical-component records of the previously identified time windows with PKP amplitude > 2 in the range ±0.05 s km −1 , using a finer slowness grid spacing of 0.0032 s km −1 .
We then convert these high-resolution PKP slowness vectors to source locations using the PKP slowness-distance relation computed with the IASP91 earth model (Fig. 3a, b; Kennett et al. (1995)).When different time windows have the same PKP slowness vector, we regard them as having the same source locations and record their cumulative duration (number of hours; Fig. 3a, b), which is sufficient for a preliminary characterization of these sources.We note that these estimated source locations are only approximate, as the slowness peaks are relatively broad in our images, 3D heterogeneity likely introduces deviations between observed slownesses and those predicted by 1D models, and the ocean-wave sources themselves are spatially defused rather than concentrated like earthquakes.A more detailed study of the spatial extent and temporal evolution of these sources will require back-projection imaging using data collected by arrays with a larger aperture than used here, which is beyond the scope of this study.
Our PKP sources are predominantly located in the Southern Ocean, where the ocean waves are the highest among all water bodies in the PKP range of our array (Fig. 3a).
We also observe far more PKP sources in the southern winter (Jul 2012-Sep 2012and Apr 2013-Sep 2013) than the southern summer (Oct 2012-Mar 2013) of our observation period (Fig. 3a), which is likely due to the greater wave height in the Southern Ocean in winter.In addition to wave height, a proxy for wave energy, P-wave radiation of oceansolid-earth interactions is also controlled by wave period and ocean depth, which can be characterized using the ocean site effect (Gualtieri et al., 2014).Our PKP sources appear to be mostly located in areas with high P-wave ocean site-effect at 4 and 5 s from Gualtieri et al. (2014) (Fig. 3b).The correlations between the spatial distribution of our PKP sources and the wave height and the ocean-site effect indicate that our PKP waves likely result from the nonlinear interaction of ocean gravity waves generated by storms (Gualtieri et al., 2014), consistent with the conclusions from previous studies that identified PKP energy in Earth's noise field (e.g.Koper andde Foy (2008), Gerstoft et al. (2008)).
We also compare the temporal variation of our PKP signals with global earthquake activity.Fig. 3c illustrates the variation of our PKP -beam amplitude over the time period of about a year, with significantly stronger PKP beams in southern winter than in southern summer.In addition to the broad peaks likely due to ocean activity, the PKPbeam amplitude shows many narrow spikes, which appear to be correlated with global seismic activity (Fig. 3c).Since the time windows containing the direct arrivals of global M > 5 events were removed from our analysis, these spikes must be due to the coda waves of the events, which can persist for hours after the first arrivals (Tkalčić et al., 2020).
Interestingly, many of these spikes correlate with events not in the PKP range (gray lines in Fig. 3c), suggesting that the coda waves of global earthquakes contain waves traveling with smaller slownesses and thus steeper incident angles than the direct phases.This observation agrees with recent studies using these steeply incident coda waves to explain the phases in Earth's correlation wave field (e.g.Tkalčić et al. (2020)).

PKS-PKP from Cross-component Cross-correlation
Wave fields dominated by a single PKP noise source are analogous to those generated by earthquakes because the wave fields in both cases are close to unidirectional.Therefore, imaging techniques designed for earthquake data, e.g., receiver function techniques, may also be applicable to noise data dominated by a single PKP noise source.
Here, we use cross-correlation between the horizontal-and vertical-component noise records as an approximation of the deconvolution procedure in receiver-function analysis (Ammon, 1991).To enhance the near-vertically traveling PKP waves while reducing surface-wave energy, which typically dominates Earth's noise field, we stack the vertical-and horizontalcomponent records of all the active stations in the array before performing cross-correlation on the stacked records (hereafter "array stacking" ).The resulting E-Z and N-Z crosscorrelation functions show a clear arrival at ∼215 s, whose amplitude appears to temporally correlate with the PKP -beam amplitude (Fig. 4a, b).This correlation is more clearly shown when we compare the temporal variation of the relative amplitude of the 215-second phase, defined as the ratio between the average absolute amplitude in a 30second window around 215 s and that in a 90-second window around 215 s, on daily stacked cross-correlation functions (red in Fig. 4b) with the temporal variation of the PKP -beam amplitude (black in Fig. 4b).This correlation suggests an association of this phase with the interaction between P waves and Earth's core.Following previous noise-imaging studies, we stacked the cross-correlation functions of many time windows to enhance the signalnoise-ratio of this phase (hereafter "215-second phase").The results show that stacking using only the time windows with a strong PKP beam produces a stronger 215-second phase than stacking using both the time windows with and without strong PKP beams (Fig. 4c-d), which is expected because the time windows without strong PKP beams generally do not show a clear 215-second phase 4a, b).Hereafter, we will focus on the time windows with PKP -beam amplitude > 2, which likely contain the highest-quality 215-second phases (Fig. 4).
To test the effects of array stacking on the waveform quality, we computed the stacked cross-correlation functions for each station individually before stacking them.Note that the difference between this method without array stacking and the method with array stacking is whether stacking across different stations is performed after (without array stacking) or before (with array stacking) cross-correlations.The comparison between the results of these two methods clearly shows that the method with array stacking produces significantly stronger 215-second phases (Fig. 4c), which is likely because stacking the noise records across the array enhances the near-vertically traveling PKP noise and the associated phases, which are responsible for the 215-second phase.From now on, we will show only the results from the cross-correlation functions with array stacking.
To further characterize the 215-second phase, we binned the PKP slowness vectors into grids with 15°and 0.005 s km −1 spacing in azimuth and slowness and stacked the cross-correlation functions of the time windows in each bin (hereafter "PKP -source bin"), which is analogous to receiver-function stacks for groups of nearby earthquakes.
While processing the data for each PKP -source bin, we aligned the records of individual stations using the back azimuth and slowness of the bin before performing array stacking, which further enhances the PKP signals.The stacked waveform shows that although the amplitude of the 215-second phase varies significantly across different source bins, its arrival time stays almost the same (Fig. 5).We also computed the best-fitting linear polarization direction for the 215-second arrival of each PKP source by finding the direction that maximizes the maximum absolute amplitude of the 215-second arrival, which is taken in a 30 s time window around 215 s, on the signal projected to the direction.These polarization directions (red bars in Fig. 5) agree very well with those of the corresponding sources bins (black bars in Fig. 5), suggesting that the 215-second phase consists of mostly SV energy.
Based on the above observations about our 215-second phase, we interpret it as PKS-PKP (Fig. 1c).Because travel-time curves of the same branches of PKP and PKS are almost parallel (Fig. 1d), the differential travel time of the two phases stays at ∼215 s across a broad range epicentral distance, which is consistent with the observation that our 215-second phase remains at approximately the same time for sources with different slownesses (Fig. 5).The radial polarization of our 215-second phase also agrees with that of PKS, which consists only of SV waves in an isotropic earth.Although different branches of PKP and PKS often arrive in the same distance range (Fig. 1c, d), the nearconstant arrival time of our 215-second phase indicates that it most likely results from the cross-correlation of PKP and PKS phases from the same branch.One possible explanation for this observation is that the different ray paths of different PKP and PKS branches leave different structural imprints on their waveform, which causes them to decorrelate.
Among our PKP beams, many have slownesses >0.032 s km −1 , which suggests that they belong to the PKPab branch.However, PKPab does not coexist with PKSab at the same distance (Fig. 1d), which appears to suggest that their clear PKS-PKP signals (e.g.Fig. 5a) result from cross-correlation between PKPab and PKS of other branches.To investigate this issue, we performed beamforming using the same dataset for four earthquakes from the USGS earthquake catalog (EQ1-4) that are close to one of our PKP sources with slowness >0.032 s km −1 (2013-07-06-00-00-00; Fig. S1).Among them, EQ1 and EQ2 show good agreement between the observed and predicted slowness, whereas EQ3 and EQ4 show greater slownesses than the predictions (Fig. S1c), which are probably due to lateral heterogeneity along the ray paths.We thus infer that our PKP beams with >0.032 s km −1 may actually represent PKPbc waves whose slownesses are elevated due to similar 3D structural effects, which, unlike PKPab, coexist with PKSbc at the same distance.We note that the 3D structural effects likely also cause errors in our PKP source locations, which should only be regarded as preliminary estimates.

Discussion
To our knowledge, this is the first report of PKS-PKP retrieved from noise data.
Although our PKS-PKP observation has the same arrival time (∼215 s) as cS-cP, a phase in Earth's correlation wavefield, at zero station offset (Pha . m et al., 2018), the two phases are fundamentally different for two main reasons: First, our PKS-PKP has its counterpart in earthquake records PKS, whereas cS-cP is not observed in earthquake records.
Second, our PKS-PKP is extracted via cross-correlation of different data components recorded at the same location, whereas cS-cP is retrieved through cross-correlation of vertical-component data recorded at different locations (Pha . m et al., 2018).Because PKS is routinely used for shear-wave-splitting analyses (e.g.Long and Silver (2009)), we also experimented with shear-wave splitting analysis (see Supplementary Text 1 for the method) using our PKS-PKP observations but obtained results very different from previous studies.The two PKP -source bins with the clearest PKS-PKP waveforms, PKP01 and PKP05 (Fig. 5), yielded fast directions of 121°and 127°, respectively (Figs.S2 and S3), significantly different from ∼70°given by shear-wave-splitting analyses of earthquake data (Yang et al., 2017).This discrepancy could be due to the low quality of our signals as the eigenvalueratio distributions indicate that neither of the two measurements is very conclusive (Figs.S2c and S3c).Since our data contain energy only in the narrow band between 2 and 10 s, whereas earthquake data typically contain more long-period energy, another possible explanation for this discrepancy is that our results are affected more by shallow structure than those from earthquake data.This hypothesis is supported by previous studies showing increased sensitivity of SKS splitting parameters to shallow structure at shorter periods (e.g.Sieminski et al. (2008)).In addition, Wirth and Long (2014) gave a NW-SE fast direction in the upper lithosphere of our study area, which is more consistent with our results.
Although the arrival time of our PKS-PKP observations stay mostly the same for different PKP sources, its amplitude varies significantly (Fig. 5).This variation does not appear to be due to stacking fold because sources with lower stacking fold can have stronger PKS-PKP than those with higher stacking folder (e.g.PKP05 compared with PKP04).
Therefore, the variation is likely due to differences in the sources or the structures that the waves travel through.The sources with stronger PKS-PKP may radiate stronger PKP waves.Alternatively, heterogeneity at the core-mantle boundary (CMB), e.g. the Ultra Low Velocity Zones (Garnero et al., 1998), may cause changes in PKS waveforms.One way to separate contributions from source and structure is to observe PKS-PKP across a broader range.The Transportable Array (TA) is suitable for this purpose, although its station density is significantly lower than that used here.Nonetheless, we may be able to achieve a similar signal quality with the TA data by stacking stations within a broader radius (the current limit is a 100 km-radius circle) because the increased range will still be much smaller than the depth to the CMB.
Our results show that PKP noise sources are extremely variable in both space and time, which likely also applies to other body-wave noise sources.We also find that bodywave scattering signals extracted from noise data can be significantly enhanced with simple techniques, namely time-window selection and array stacking, that address the spatiatemporal variation of body-wave sources.In principle, time-window selection does not require dense-array data, although a synchronous array may be necessary to determine the time windows containing significant body-wave noise energy.Array stacking requires array data, which limits its application, although the required array density likely depends on the targeted seismic phase.So far, most of the seismic imaging studies using body-wave noise have not accounted for its spatiatemporal variation and have relied simply on stacking large number of cross-correlation functions (e.g.Poli et al. (2012) and Feng et al. (2021)).Our results suggest that the primary contribution to their signals may have only come from a fraction of all the time windows, and that simply selecting those time windows might significantly improve the signal quality (Fig. 4).The signal quality may be further improved if array stacking can be performed before cross-correlation.

Conclusions
We extract a phase that likely represents PKS-PKP from cross-component crosscorrelation of noise recordings.We show that the amplitude of PKS-PKP is significantly enhanced when only time windows containing strong PKP signals are used.We also show that stacking array data before cross-correlation significantly enhances PKS-PKP amplitudes.Future studies that retrieve body-wave scattered phases from noise data should account for the spatiatemporal variation of body-wave noise sources.

Data Availability Statement
The seismic and wave-height data used in this study are freely available through the Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC) https://ds.iris.edu/ds/nodes/dmc/and the Environmental Modeling Center of NOAA https://polar.ncep.noaa.gov/waves/wavewatch/,respectively.The plots in this paper are created with the Generic Mapping Tools (Wessel et al., 2019).

Acknowledgments
This study is funded by NSF Grants EAR-1358510 and EAR-1829601.T.L. is supported by a Green Postdoctoral Scholarship.IRIS DMC is funded by the the NSF under Cooperative Support Agreement EAR-1851048.We thank Lucia Gualtieri for providing the ocean site-effect maps and Wenyuan Fan for stimulating discussion.
Most of these studies focused on extracting surface-wave signals from the noise field because surface waves usually dominate the signals retrieved by noise cross-correlation.This observation is commonly attributed to the prominence of surface waves in Earth's noise field as a result of noise sources, such as wind and ocean waves, occurring mostly at the surface.Despite their lower amplitudes, body-wave signals have occasionally been retrieved from noise cross-correlations and used to image Earth structure (e.g., Feng et al., 2021;Nakata et al., 2015;Pedersen and Colombi, 2018;Poli et al., 2012).A major advantage of body waves over surface waves in studying Earth structure is that body-wave reflected and converted phases are sensitive to material discontinuities in Earth's interior (e.g., the Moho and the core-mantle boundary [CMB]), which cannot be resolved with surface-wave data alone.However, body-wave reflection and conversion signals are weaker than direct phases and thus more difficult to observe in the cross-correlation functions, which are typically nosier than earthquake records.Therefore, techniques capable of enhancing body-wave reflection and conversion signals are needed to better image Earth's discontinuities with noise records.
In addition to imaging using seismic noise, in recent years, major advances have been made in understanding the sources of Earth's noise field (e.g., Gualtieri et al., 2014;Liu et al., 2020;Nishida and Takagi, 2016;Retailleau and Gualtieri, 2021).Many contributions were made by studying body-wave noise signals with array techniques (e.g., beamforming and back-projection), which suggests that weak body-wave noise signals can be enhanced with array processing to better image Earth structure.These studies also showed that body-wave noise sources, which are usually associated with storms in the oceans, are likely spatially and temporally heterogeneous, which Abstract Seismic noise has been widely used to image Earth's structure in the past decades as a powerful supplement to earthquake signals.Although the seismic noise field contains both surface-wave and body-wave components, most previous studies have focused on surface waves due to their large amplitudes.Here, we use array analyses to identify body-wave noise traveling as PKP waves.We find that by cross-correlating the array-stacked horizontal-and vertical-component data in the time windows containing the PKP noise signals, we extract a phase likely representing PKS-PKP, the differential phase between PKS and PKP.This phase can potentially be used for shear-wave-splitting analysis.Our results also suggest that the sources of body-wave noise are extremely heterogeneous in both space and time, which should be accounted for in future studies using body-wave noise to image Earth structure.
Plain Language Summary Seismic noise is the vibration of Earth generated by activities other than earthquakes, such as wind and ocean waves.Signals extracted from seismic noise can be used to study Earth's interior structure in ways similar to how earthquake records have been analyzed.Most previous studies using seismic noise to study Earth structure used its surface-wave component, that is, the waves propagating at Earth's surface, whereas the body-wave component, that is, the waves traveling through Earth's interior, is less used because body-wave noise is usually much weaker than surface-wave noise.Here, we use data collected by a dense seismic array to identify body-wave noise propagating as PKP waves, P waves that travel through Earth's core.We also find that PKS-PKP, the differential phase between PKS and PKP, can be extracted from the records of time windows containing strong PKP energy.This phase can potentially be used to study the anisotropic properties of Earth's crust and mantle.Here, we present observations of body-wave noise propagating as PKP using data collected by a dense broadband seismic array in the central US.We further show that a phase likely representing PKS-PKP can be extracted by cross-correlating the array-stacked horizontal-and vertical-component noise records in the time windows containing the PKP noise signals.We then discuss the potential applications of this phase and the implications of our findings for seismic interferometry.

Data and Preprocessing
We mainly use the continuous data collected by the Ozark Illinois Indiana Kentucky (OIINK) Flexible Array Experiment (network code: XO), a dense 2D broadband seismic array with a station spacing of ∼25 km located in the central US (Figures 1a and 1b).To make the resolution of our results more isotropic, we select the OIINK stations located in a 100-km radius circle and also include the Transportable-Array stations in this range (Figure 1b).Although the two arrays together span 2011-2015, to ensure a reasonable resolution, we focus on time windows with more than 20 active stations, which limits our analysis to a roughly 1-year period between June 2012 and August 2013.We downloaded the continuous data from the IRIS Data Management Center in 1-hr time windows, removed the instrument response, and band-pass filtered the data to 2-10 s, which contains the secondary microseism energy.To avoid the effects of earthquakes and instrument malfunctions, we removed the 1-hr time windows containing the first arrivals of global earthquakes with magnitude >5 and those containing amplitudes >1 × 10 −15 m s −1 .(Kennett et al., 1995).Inset shows the differential travel times between PKS and PKP of the same branches.

PKP Signals From Beamforming Analysis
We performed conventional linear beamforming with all three components (vertical, east, and north) of our array data to characterize the directional properties of the noise field.To save computational cost, we first performed a reconnaissance analysis over the slowness range of ±0.2 s km −1 at a grid spacing of 0.013 s km −1 in the W-E and S-N directions.The resulting vertical-component slowness images clearly show beams with slowness <0.04 s km −1 (Figure 2a), which likely represent PKP signals as suggested by previous studies (Koper & de Foy, 2008;Landès et al., 2010).The horizontal-component slowness images also show local maxima corresponding to the PKP beams on the vertical component, though the background noise is significantly higher on the horizontal-component images (Figure 2a), which could be due to the near-vertical particle motion of PKP or a more homogeneous distribution of horizontal-component noise sources.The slowness images of some time windows also show multiple peaks (e.g., 2013-07-06-00-00-00; Figure 2a).
For seismic imaging, we prefer to use time windows dominated by PKP energy from a single direction because this source distribution resembles that of earthquake sources, which may make techniques in earthquake imaging readily applicable.To identify these time windows, we find the maximum in the PKP range (slowness <0.04 s km −1 ) of each vertical-component slowness image and the corresponding slowness vector, which we refer to as the PKP slowness.We then define the vertical-component normalized PKP-beam amplitude as the ratio between the maximum amplitude in the PKP range and the average amplitude of the whole slowness image, which measures the power of the strongest PKP beam relative to the background noise.We further define the corresponding normalized PKP-beam amplitudes for the horizontal components as the ratios between the amplitudes at the PKP slowness measured previously from the vertical-component slowness image and the average amplitudes of the whole slowness images.We finally define the three-component normalized PKP-beam amplitude (hereafter "PKP-beam amplitude") as the product of the normalized PKP-beam amplitudes of the three components.We regard the time windows with PKP-beam amplitude >2, which account for about 10% of all the time windows, as windows dominated by PKP energy from a single direction (hereafter "PKP windows").
To enhance the slowness and back-azimuth resolution for the PKP beams, we further performed beamforming for the vertical-component records of the PKP windows in the range ±0.05 s km −1 , using a finer grid spacing of 0.0032 s km −1 .A histogram of the resulting slownesses shows that the vast majority of these time windows are dominated by PKPbc beams close to the b caustic (Figure 2b), which is likely due to the amplification of PKP near its caustics.A significant number of windows show slownesses >0.032 s km −1 , which suggests that they are dominated by PKPab beams.However, beamforming results of earthquakes with known locations near these sources indicate that these apparent PKPab beams are probably PKPbc beams with elevated slownesses due to the effects of 3D velocity structure (Supplementary Text 1 and Figure S1 in Supporting Information S1).
Our PKP-beam amplitude shows a clear seasonal variation, with high amplitude in southern winter (April-October) and low amplitude in southern summer (November-March; Figure 2c).This seasonality is likely due to higher waves in the Southern Ocean in southern winter, where most of the PKP energy is generated through ocean-solid-earth interaction (Supplementary Text 2 and Figure S2 in Supporting Information S1).Interestingly, our PKP-beam amplitude also shows some narrow spikes that correlate with global earthquake activities (Figure 2c).Since the time windows containing the direct arrivals of global M > 5 events were removed from our analysis, these spikes must be due to the late coda waves of these events, which can persist for hours after the first arrivals (Tkalčić et al., 2020).Many of these spikes correlate with events not in the PKP range (gray lines in Figure 2c), suggesting that the coda waves of global earthquakes contain waves traveling with smaller slownesses and thus steeper incident angles than the direct phases.This observation agrees with recent studies using these steeply incident coda waves to explain the phases in Earth's correlation wavefield (e.g., Tkalčić et al., 2020).We also find the source locations for PKP beams outside the late-coda windows (Supplementary Text 2 in Supporting Information S1), which agree well with the significant wave-height data from WAVEWATCH III (Tolman, 2009) (Figure S2a in Supporting Information S1) and the ocean site effect map from Gualtieri et al. (2014) (Figure S2b in Supporting Information S1).

PKS-PKP From Cross-Component Cross-Correlation
Wave fields dominated by a single PKP noise source are analogous to those generated by earthquakes because the wave fields in both cases are close to unidirectional.Therefore, imaging techniques designed for earthquake data, for example, receiver-function techniques, may also be applicable to records in our PKP windows.Here, we use cross-correlation between the vertical-and horizontal-component noise records as an approximation of the deconvolution procedure in receiver-function analysis (Ammon, 1991).To enhance the near-vertically traveling PKP waves while reducing surface-wave energy, which typically dominates Earth's noise field, we stack the vertical-and horizontal-component records of all the active stations in the array before performing cross-correlation on the stacked records (hereafter "array stacking").Our initial stacking of the entire data set assumes zero slowness (i.e., vertical wave propagation); later we will refine our stacks to sum more accurately the energy seen arriving at particular back azimuths and slownesses during specific time intervals.Note that stacking before and after performing cross-correlation are different because the former also includes the cross terms between different stations (Supplementary Text 3 in Supporting Information S1).
Our E-Z and N-Z cross-correlation functions show a clear arrival at ∼215 s, whose amplitude appears to temporally correlate with the PKP-beam amplitude (Figures 3a and 3b).This correlation is more clearly shown when we compare the temporal variation of the relative amplitude of the 215-s phase, defined as the ratio between the average absolute amplitude in a 30-s window around 215 s and that in a 90-s window around 215 s, on daily stacked cross-correlation functions (red in Figure 3b) with the temporal variation of our PKP-beam amplitude (black in Figure 3b).This correlation suggests an association of this phase with the interaction between P waves and Earth's core.Following previous noise-imaging studies, we stacked the cross-correlation functions of many time windows to enhance the signal-noise ratio of this phase (hereafter "215-s phase").The results show that stacking using only the time windows with a strong PKP beam produces a stronger 215-s phase than stacking using both the time windows with and without strong PKP beams (Figures 3c-3e), which is expected because the time windows without strong PKP beams generally do not show a clear 215-s phase (Figures 3a and 3b).Hereafter, we will focus on our PKP windows (time windows with PKP-beam amplitude >2), which likely contain the highest-quality 215-s phases (Figure 3c).
To test the effects of array stacking on the waveform quality, we also compared the results with and without array stacking, which clearly shows that the method with array stacking produces significantly stronger 215-s phases (Figure 3c).This is likely because stacking the noise records across the array enhances the near-vertically traveling PKP noise and its associated phases, which are responsible for the 215-s phase.From now on, we will show only the results with array stacking.
Since the time windows that we used to extract the 215-s phase also include windows containing global-earthquake late coda (<10 hr after the events; Figure 2c), an important question is whether the main contribution of our 215-s phase comes from earthquake coda energy.To investigate this possibility, we excluded the time windows <10 hr after global M > 5 events and performed the same analysis.The results show that despite the removal of nearly 3/4 of the original time windows, the 215-s phase remains clear on the stacked cross-correlation function, though with slightly lower signal-noise ratio due to the lower stacking fold (Figure S3 in Supporting Information S1).Moreover, the stack including only the time windows with strong PKP beams still shows a stronger 215-s phase than the one including all time windows (Figures S3b-S3d in Supporting Information S1).These results clearly demonstrate that global-earthquake late coda is not the only cause of our 215-s phase, with ocean-solidearth interaction likely also contributing significantly as evidenced by the clear seasonality of our PKP-beam amplitude (Figure 2c).We note that our data have a period band (2-10 s) much shorter than data typically used for earthquake-late-coda analyses (>15 s; e.g., Wang and Tkalčić, 2020).Boué et al. (2014) demonstrated that in our short-period band, noise cross-correlations are largely unaffected by earthquake late coda, probably because the coda waves generally lack short-period components due to the high cumulative attenuation along their long paths, although some events may be more efficient in generating short-period signals, which cause the PKP-energy bursts that correlate with global seismic activities (Figure 2c).We thus conclude that global earthquake coda does not contribute significantly to our 215-s phase.
To further characterize our 215-s phase, we binned the slowness vectors of our PKP windows into grids with 15° and 0.005 s km −1 spacing in azimuth and slowness, respectively (hereafter "PKP-source bin"; Figure 4).Since the PKP waves of these source bins have small yet nonzero slownesses (first column of Figure 4), stacking the noise records without applying time shifts, which is equivalent to assuming zero slowness, will not maximize the energy of PKP and its secondary phases, which appear correlated with our 215-s phase (Figure 3b).To find the slowness vectors that maximize the amplitude of our 215-s phase, we performed array stacking assuming a range of slowness vectors for each source bin and found the 215-s-phase amplitude on the stacked cross-correlation functions for each slowness vector, which is defined as the maximum amplitude in the time window 200-240 s.This procedure gives a 215-s-phase-amplitude slowness image (hereafter "phase-amplitude image") for each source bin (second column of Figure 4).
The phase-amplitude images of PKP02, 04, and 05 clearly show maxima at slowness vectors very similar to the corresponding beam-amplitude images (Figures 4b,4d,and 4e), which indicates that performing array stacking after shifting the noise records by the PKP slowness enhances the 215-s phase more than stacking without time shifts.The phase-amplitude image of PKP01 shows a diffuse energy distribution, and the phase-amplitude image of PKP03 shows a maximum at a significantly smaller slowness than the beam-amplitude image (Figures 4a and 4c).We will later show that these discrepancies are likely due to the unusually large slownesses of the two bins.To maximize the amplitude of our 215-s phase, we thus shifted the noise records of PKP01-05 using their PKP slownesses before array stacking and cross-correlation.The stacked cross-correlation functions of PKP02, 04, and 05 show clear 215-s phases, whereas the phase is less visible on PKP01 and 03 (third column of Figure 4).We also computed the best-fitting linear particle-motion direction for the 215-s phase of each PKP-source bin.For PKP01, 02, and 04, the linear particle-motion directions (red bars in the fourth column of Figure 4) agree well with the back azimuths of the corresponding source bins (black bars in the fourth column of Figure 4), suggesting that the 215-s phase consists of mostly SV energy.For PKP03 and 05, the linear particle-motion directions are significantly different from the source directions.We will later show that these differences are likely due to shear-wave splitting.
Based on the above observations of the 215-s phrase, we interpret it as PKS-PKP, the differential phase between PKS, the core phase with a P-to-S conversion at the receiver-side CMB, and PKP (Figure 1c).The IASP91-predicted differential travel time between PKSdf and PKPdf is 212 s at ∼180° (i.e., vertical incidence), whereas the differential travel time between PKSbc and PKPbc is ∼215 s in 140°-150° (Figure 1d).The 215-s-phase arrival times of PKP02, 04, and 05, the three source bins with high 215-s-phase amplitude, are all slightly larger than 212 s (gray vertical lines in the third column of Figure 4), which suggests that the PKP and PKS rays likely belong to the bc branch and thus are obliquely incident, consistent with the beam-amplitude and phase-amplitude images (first and second columns of Figure 4).This interpretation is supported by the fact that P-to-S conversions are not predicted for vertically traveling P waves for 1D Earth models.Note that the near-radial polarization of the 215-s phase of PKP01, 02, and 04 also agrees with that of PKS (fourth column of Figure 4), which consists only of SV waves in an isotropic earth.Hereafter, we will use PKP-PKS to refer to our 215-s phase.PKP01 and 03 have slownesses greater than the b caustic of PKP given by IASP91 probably due to the effects of 3D Earth structure (Figures 2b, 4a and 4c).Because the slowness difference between PKPbc and PKSbc recorded at the same distance grows with increasing slowness (Figure 2b), the unusually large slownesses of PKP01 and 03 may cause their PKP and PKS to have sufficiently different slownesses that the two phases cannot be enhanced with array stacking using one single slowness vector, which could explain the low PKS-PKP amplitude and the discrepancy between the beam-amplitude and phase-amplitude images of PKP01 and 03 (Figures 4a and 4c) In theory, SKP, the core phase with S-to-P conversion at the source-side CMB, arrives at the same time as PKS and thus might cause interference, assuming S waves are generated at the source region.However, since SKP arrives as a P wave at the receiver, it likely has very low amplitude on the horizontal components due to its near-vertical rays, which should make it much weaker than PKS on our vertical-horizontal cross-correlation functions.
Because PKS is routinely used for shear-wave-splitting studies (e.g., Long and Silver, 2009), we also performed shear-wave-splitting analysis (see Supplementary Text 4 in Supporting Information S1 for the method) on our PKS-PKP observations.Among PKP01-05, PKP05 yields the best shear-wave-splitting results as evidenced by its diagnostic elliptical particle motion before time correction (Figure 5b) and well-focused maximum on the eigenvalue-ratio distribution (Figure 5c).The fast-direction (46°) and splitting time (1.4 s) are reasonably close to the results of Yang et al. (2017) derived from earthquake data recorded at stations located within our circular array window (Figure 5c).We also derive a similar set of splitting parameters (66° and 1.4 s) from PKP03 (Figure S4 in Supporting Information S1), a bin with a source direction nearly orthogonal to that of PKP05 (Figure 4c), though the splitting parameters are less well constrained likely due to the low amplitude of PKS-PKP.The other source bins produce only ambiguous results.Assuming a fast direction of ∼55° in our study region, PKP01, 02, and 04 all have source directions close to either the fast or the slow directions, which likely causes them to show little shear-wave splitting and near-radial particle motions (Figures 4a, 4b, and 4d).In contrast, PKP03 and 05 have source directions significantly different from both the fast and slow directions, which causes them to show significant splitting and non-linear particle motions (Figures 4c and 4e).In summary, a fast direction of ∼55° is consistent with our observations.Since our shear-wave splitting results can be regarded as derived from only two sources, whereas the ones from Yang et al. (2017) are the average results of many sources, the difference between the two might be due to lateral variation of anisotropy beneath the study region, which can cause differences between different ray paths.Another possibility is that our results are affected more by shallow structure than those from earthquake data because our data contain energy only in the short period band of 2-10 s, whereas earthquake data typically contain more long-period energy.This hypothesis is supported by previous studies showing increased sensitivity of SKS splitting parameters to shallow structure at shorter periods (e.g., Sieminski et al., 2008).These issues warrant further study, including detailed comparisons at individual stations between 10.1029/2021GL097034 9 of 11 shear-wave splitting results derived from earthquake records and those obtained from noise cross-correlation.However, in regions that PKS-PKP can be observed from noise, analysis of this phase should help contribute to upper-mantle anisotropy studies.

Discussion
To our knowledge, this is the first report of PKS-PKP retrieved from noise data.Our PKS-PKP observation can be regarded as belonging to the same broad category as the phases in the Earth's correlation wavefield (e.g., Pham et al., 2018;Tkalcic et al., 2020;Wang and Tkalčić, 2020), which are also produced through cross-correlation of records rich in steeply incident body-wave energy (global earthquake coda waves).Specifically, our PKS-PKP has an arrival time close to cS-cP at zero offset (e.g., Pham et al., 2018).Nonetheless, the generation mechanism of our PKS-PKP observations is likely different from that of cS-cP, which causes the two phases to have different structural sensitivities.cS-cP is thought to be formed via the interference between any earthquake-coda phase pairs with one P-S differential leg between the CMB and the surface (e.g., Figure 3b in Wang and Tkalčić (2020)).Because many different phase pairs satisfy this condition, and that the P-S differential leg can constitute any segment on the ray paths of the two interfering phases, the waveform of cS-cP is likely sensitive to Earth structure in a very broad range (Wang & Tkalčić, 2020).In contrast, our PKS-PKP arises mostly from PKP and PKS excited by ocean-solid-earth interactions, providing sensitivity mainly to the mantle structure beneath the station.
This interpretation is supported by four lines of evidence: First, the short period band that we focus on (2-10 s) is known to be largely free from the effects of earthquake late coda (Boué et al., 2014).Second, our PKS-PKP likely represents incoming S waves at the receiver because it is extracted via vertical-horizontal cross-component cross-correlation.Therefore, the P-S differential leg of our PKS-PKP must be the last leg on the ray paths of the two interfering phases.Third, our array stacking enhances incoming PKP and PKS waves with a specific slowness vector while suppressing contributions from phase pairs with different slowness vectors, for example, possible earthquake late coda with steeper incident angles (Figure 4).Finally, the clear shear-wave-splitting signal of PKP05, which agrees with previous results, indicates that our PKS-PKP is indeed primarily sensitive to the structure immediately below the station.In summary, our PKS-PKP is related to yet different from cS-cP and other phases in Earth's correlation wavefield.
Our results show that PKP noise sources are extremely variable in both space and time, which likely also applies to other body-wave noise sources.We also find that body-wave scattering signals extracted from noise data can be significantly enhanced with simple techniques, namely time-window selection and array stacking, that address the spatiotemporal variation of body-wave sources.In principle, time-window selection does not require dense-array data, although a synchronous array may be necessary to determine the time windows containing significant body-wave noise energy.Array stacking requires array data, which limits its application, although the required array density likely depends on the targeted seismic phase.So far, most of the seismic imaging studies using body-wave noise have not accounted for its spatiatemporal variation and have relied simply on stacking large number of cross-correlation functions (e.g., Feng et al., 2021;Poli et al., 2012).Our results suggest that the primary contribution to their signals may have only come from a fraction of all the time windows, and that simply selecting those time windows might significantly improve the signal quality (Figure 3).The signal quality may be further improved if array stacking can be performed before cross-correlation.

Conclusions
We extract a phase that likely represents PKS-PKP from cross-component cross-correlation of noise records.We show that the amplitude of PKS-PKP is significantly enhanced when only time windows containing strong PKP signals are used.We also show that stacking array data before cross-correlation significantly enhances PKS-PKP amplitudes.The shear-wave-splitting parameters estimated with our PKS-PKP waveforms are similar to the ones from previous studies derived with earthquake data, suggesting that PKS-PKP may be used for studying crust and mantle anisotropy in the future.

Data Availability Statement
The metadata of TA and XO can be accessed at https://ds.iris.edu/mda/TA/and https://ds.iris.edu/mda/XO/?start-time=2011-01-01T00:00:00&endtime=2015-12-31T23:59:59,respectively.The time-series data of the two arrays are freely availabe at the Incorporated Research Institutions for Seismology Data Management Center and were downloaded using ObsPy in this study (Krischer et al., 2015).The wave-height data of WAVEWATCH III are freely available at the Environmental Modeling Center of NOAA (https://polar.ncep.noaa.gov/waves/wavewatch/).The P-wave site-effect maps in this paper are provided by Lucia Gualtieri through personal communication and are available at https://doi.org/10.5281/zenodo.5904118.The plots in this paper are created with the Generic Mapping Tools (Wessel et al., 2019).

Figure 1 .
Figure 1.Station locations and PKP and PKS ray geometries and travel times.(a) Map of the contiguous US showing the closeup of panel (b) marked in red.(b) Map of all the OIINK stations (magenta) and nearby TA stations (cyan).The 100-km radius circle defines the region in which the stations are included in our analysis.(c) Ray paths of PKPab, PKPdf (blue), and PKSdf (cyan) at 160°.(d) Travel times as functions of epicentral distance for different branches of PKP (blue) and PKS (cyan)

Figure 3 .
Figure3.Spatial distribution and temporal variation of our PKP sources.(a) Spatial distributions of our PKP sources overlain on the ocean site-effect maps for period = 4 s (left) and 5 s (right) fromGualtieri et al. (2014).Sizes of the circles denote the cumulative duration of each source.(b) The same as (a), but for sources in the southern winter (left) and summer (right) of our observation period overlain on the average significant wave-height maps for the respective seasons from WAVEWATCH III(Tolman et al., 2009).(c) Three-component PKP -beam amplitude as a function of time.Red and gray lines mark the origin times of global M > 6 events in and out of the PKP epicentral distance range, respectively.

Figure 4 .
Figure 4. E-Z and N-Z cross-correlation functions of the array-stacked records.(a) E-Z (left) and N-Z (right) cross-correlation functions for all the active time windows in a three-month period from June to September 2012.(b) Temporal variation of PKP -beam amplitude (black) and the 215-second-phase amplitude (red) for the time range in (a).(c-d) Blue waveform: Stacked E-Z (left) and N-Z (right) cross-correlation functions for time windows with PKP -beam amplitude (c) > 2, (d) > 1, and (e) all the time windows.Gray waveform in (c): The same as the blue waveform, but computed with stacking E-Z and N-Z cross-correlation functions of individual stations.

Figure 5 .
Figure 5. Stacked cross-correlation functions for the five PKP-source bins with the most cumulative duration: (a) PKP01, (b) PKP02, (c) PKP03, (d) PKP04, and (e) PKP05.Left column: Stacked E-Z (blue) and N-Z (yellow) cross-correlation functions.Right column: PKP beam direction (black) and the best-fit linear polarization (red) for the signals in a 30-s time window around 215 s.

Figure S1 .
Figure S1.Using earthquakes with known locations to evaluate biases in our beamforming.(a) Locations of Earthquake (EQ) 1-4 used for calibration (white stars) and the derived PKP-source location of time window 2013-07-06-00-00-00 (yellow star).(b) Slowness image of the time window 2013-07-06-00-00-00 with the maximum marked with a dark blue cross.The gray circle denotes the slowness of 0.04 s km -1 .(c) Same as (b), but for EQ 1-4.The light blue crosses mark the slowness vectors predicted with IASP91 (multiple slownesses are due to different PKP branches).

Figure S2 .Figure S3 .
Figure S2.PKP source locations of our PKP windows that are less likely affected by earthquake late coda.Circle size denotes the cumulative duration of a certain location.(a) Sources in southern winter (left) and southern summer (right) plotted on the average significant waveheight maps of the corresponding seasons from WAVEWATCH III.(b) Sources plotted on the ocean P-wave site-effect maps at 4-s and 5-s periods from Gualtieri, et al., 2014.

Figure S4 .
Figure S4.Same as Fig. 5, but for the source bin PKP03.

Figure 1 .
Figure 1.Station locations and PKP and PKS ray geometries and travel times.(a) Map of the contiguous US showing the closeup of panel (b) marked in red.(b) Map of all the Ozark Illinois Indiana Kentucky stations (magenta) and nearby TA stations (cyan).The 100-km radius circle defines the region in which the stations are included in our analysis.(c) Ray paths of different travel-time branches of PKP (blue) and PKS (cyan) at 145°.(d) Travel times as functions of epicentral distance for different branches of PKP (blue) and PKS (cyan) computed using IASP91(Kennett et al., 1995).Inset shows the differential travel times between PKS and PKP of the same branches.

Figure 2 .
Figure 2. PKP beams characterized with array analyses.(a) Example three-component slowness images for two one-hour time windows 2013-07-06-00-00-00 (top) and 2012-07-16-11-00-00 (bottom) with clear PKP energy.Gray circle: slowness of 0.04 s km −1 .(b) Slowness-distance relation of PKP (blue curve) and PKS (cyan curve), and the slowness histogram of the PKP windows.(c) Three-component PKP-beam amplitude as a function of time.Red and gray lines mark the origin times of global M > 6 events in and out of the PKP epicentral-distance range, respectively.

Figure 3 .
Figure 3. E-Z and N-Z cross-correlation functions.(a) E-Z (left) and N-Z (right) cross-correlation functions computed with array stacking (no time shifts applied to individual traces) for all the active time windows in a 3-month period from June to September 2012.(b) Temporal variation of PKP-beam amplitude (black) and the 215-s-phase amplitude (red) for the time range in (a).(c-d) Blue waveform: Stacked E-Z (left) and N-Z (right) cross-correlation functions computed with array stacking for time windows with PKP-beam amplitude (c) >2, (d) >1.5, and (e) all time windows.Gray waveform in (c): The same as the blue waveform, but computed without array stacking.

Figure 4 .
Figure 4. Characterization of the noise recordings and cross-correlation functions of the five PKP-source bins with the longest cumulative durations: (a) PKP01, (b) PKP02, (c) PKP03, (d) PKP04, and (e) PKP05 (highlighted due to its best shear-wave splitting results).First column: Stacked beam-amplitude slowness images.Gray circles: slowness of 0.04 s km −1 .Second column: Maximum amplitudes of the 215-s phase as functions of slowness vectors used for array stacking.Third column: Stacked E-Z (blue) and N-Z (yellow) cross-correlation functions computed with array stacking using the slowness vectors of the source bins.Gray vertical line marks the arrival time of PKS-PKP at ∼180° (i.e., vertical incidence).Fourth column: Back azimuths of the source bins (black) and the best-fit linear particle-motion directions (red) of the 215-s phases.

Figure 5 .
Figure 5. Shear-wave splitting results for the PKP-source bin PKP05.(a) E-Z and N-Z components of the PKS-PKP phase before (top) and after (bottom) the time correction.(b) Particle-motion diagrams of the PKS-PKP phase before (top) and after (bottom) the time correction.(c) Normalized eigenvalue-ratio of the particle motions after time correction computed using various fast directions and splitting times.Blue cross marks the maximum.Cyan crosses mark the splitting parameters of individual stations located in our circular array window from Yang et al. (2017).