Exact free oscillation spectra, splitting functions and the resolvability of Earth's density structure

Exact free oscillation spectra, splitting functions and the resolvability of Earth's density... Summary Seismic free oscillations, or normal modes, provide a convenient tool to calculate low-frequency seismograms in heterogeneous Earth models. A procedure called ‘full mode coupling’ allows the seismic response of the Earth to be computed. However, in order to be theoretically exact, such calculations must involve an infinite set of modes. In practice, only a finite subset of modes can be used, introducing an error into the seismograms. By systematically increasing the number of modes beyond the highest frequency of interest in the seismograms, we investigate the convergence of full-coupling calculations. As a rule-of-thumb, it is necessary to couple modes 1–2 mHz above the highest frequency of interest, although results depend upon the details of the Earth model. This is significantly higher than has previously been assumed. Observations of free oscillations also provide important constraints on the heterogeneous structure of the Earth. Historically, this inference problem has been addressed by the measurement and interpretation of splitting functions. These can be seen as secondary data extracted from low frequency seismograms. The measurement step necessitates the calculation of synthetic seismograms, but current implementations rely on approximations referred to as self- or group-coupling and do not use fully accurate seismograms. We therefore also investigate whether a systematic error might be present in currently published splitting functions. We find no evidence for any systematic bias, but published uncertainties must be doubled to properly account for the errors due to theoretical omissions and regularization in the measurement process. Correspondingly, uncertainties in results derived from splitting functions must also be increased. As is well known, density has only a weak signal in low-frequency seismograms. Our results suggest this signal is of similar scale to the true uncertainties associated with currently published splitting functions. Thus, it seems that great care must be taken in any attempt to robustly infer details of Earth's density structure using current splitting functions. Computational seismology, Seismic tomography, Surface waves and free oscillations, Theoretical seismology 1 INTRODUCTION Our understanding of the Earth's large-scale interior structure and dynamics draws heavily on observations of seismic free oscillations (‘normal modes’). In particular, measurements of free-oscillation spectra at long periods represent one of the few classes of geophysical observable with meaningful sensitivity to lateral variations in density within the Earth. As such, free oscillation spectra have played a central role in attempts to determine the nature of the ‘Large Low Shear Velocity Provinces’ that have been identified within the lowermost mantle, and in attempts to identify the relative importance of thermal and chemical effects as driving forces for mantle convection (e.g. Ishii & Tromp 1999; Trampert et al.2004; Lay 2007; Mosca et al.2012; Moulik & Ekström 2016; Koelemeijer et al.2017). However, as this paper will show, computational limitations have led most such studies to rely upon seismological approximations that might be inadequate for imaging density variations. As such, their conclusions must be approached with considerable caution. In order to derive information about Earth's interior structure from seismic data, we must search for models which yield synthetic (predicted) observables which agree with those actually measured. Many strategies exist for performing this search, but all rely on the presumption that synthetic observables are computed accurately. The theoretical basis for computing synthetic normal mode spectra is well-developed, relying on a technique known as ‘normal mode coupling theory’ (e.g. Dahlen 1968, 1969; Woodhouse & Dahlen 1978; Woodhouse 1980; Woodhouse & Girnius 1982; Park 1986, 1987, 1990; Romanowicz 1987; Tromp & Dahlen 1990; Lognonné 1991; Hara et al.1991, 1993; Um & Dahlen 1992; Deuss & Woodhouse 2001; Al-Attar et al.2012; Yang & Tromp 2015; Lau et al.2015). Unlike some earlier studies, our calculations are based on the work of Al-Attar et al. (2012) and take the full effect of 3-D density variations into account, following the theory developed in Woodhouse (1980). Our calculations also allow us to take attenuation fully into account, although in the examples shown in this paper, only the spherically symmetric attenuation profile of PREM (Dziewonski & Anderson 1981) is used. Thus, our framework is physically complete, and allows the seismic response of any Earth-like body to be expressed exactly. However, the resulting expressions involve an infinite series expansion. This must be truncated in any computational implementation of the theory, introducing a truncation error within the synthetic spectra. Various common truncation schemes exist, motivated from physical principles. Computational costs scale unpleasantly with the number of terms retained in the series, so early studies had little choice but to adopt simplifying strategies. The most limiting approximation is known as ‘self-coupling’, while ‘group-coupling’ incorporates additional terms (we explain the distinction in more detail below). However, a number of studies, starting with Deuss & Woodhouse (2001), have demonstrated that neither can be relied upon within the frequency range of interest. In particular, Al-Attar et al. (2012) highlighted that use of either the self- or group-coupling approximations when computing synthetic spectra leads to errors of similar magnitude to the signal likely to be attributable to lateral variations in density within the Earth. Subsequently, a comprehensive study by Yang & Tromp (2015) concluded that self-coupling is ‘marginally acceptable’ when considering spectra below around 1.5 mHz, and ‘unacceptable’ above this point. Group-coupling was found to also yield significant errors, albeit at a lower level than in self-coupling. In principle, one could envisage improving group-coupling strategies by computing coupling strengths between modes—perhaps based on the work of Park (1987)—and ensuring that all significant interactions are accounted for. However, the notion of coupling strength is itself based on single-scattering approximations, and it is not clear how accurate the results would be. The optimal approach is therefore to adopt a strategy known as ‘full-coupling’. Despite the name, this continues to involve a truncation of the infinite series, but whereas self- and group-coupling are based around the presumption that the series is dominated by only a few terms, full-coupling aims to retain all terms that have any potential to contribute significantly; again, a more precise definition will be given in due course. Building on the existing body of work, the first aim of the present paper is to identify more precisely the conditions under which full-coupling can be regarded as ‘sufficiently accurate’. This information then permits us to perform high-quality synthetic experiments to address the second question: can robust information about earth structure be inferred from measurements of ‘splitting functions’ (Woodhouse & Giardini 1985; Giardini et al.1987, 1988)? These are a particular representation of information derived from seismic spectra, and they are an important ingredient in many current models of Earth structure (e.g. Ritsema et al.1999; Ishii & Tromp 1999; Trampert et al.2004; Ritsema et al.2011; Mosca et al.2012; Moulik & Ekström 2016; Koelemeijer et al.2017). However, the theory that connects splitting functions to structural parameters relies upon precisely the same assumptions as are inherent to self- and group-coupling. Given that these have been demonstrated to introduce significant bias into synthetic spectra, are splitting function inferences meaningful? The answer is not a foregone conclusion—the validity of any approximation is entirely context-specific—but the question must clearly be addressed. The first to address this question were Resovsky & Ritzwoller (1998). Being limited in computational resources, they estimated the theoretical error by a transfer function. Having at our disposal an arbitrarily precise forward solver, we can, for the first time, quantify this theoretical error properly. We therefore perform a simple, but definitive test: we compute accurate synthetic seismograms for the earth model S20RTS (Ritsema et al.1999), making use of full-coupling. We then measure splitting functions from this data set as if it were observed data. We also compute the ‘true’ splitting functions, by calculating the appropriate radial integrals for the known input model. Comparing the two, we find noticeable discrepancies. Repeating the experiment for synthetic spectra using the self-coupling approximation allows us to separate the contribution from theoretical approximations and regularization in the uncertainty of the measured splitting functions. Finally repeating the experiment for models where density perturbations are set to zero, we conclude that splitting functions are contaminated by a significant approximation error, a similar regularization error and contain a density signal of the level of the total measurement uncertainty. 2 THEORETICAL BACKGROUND Imaging the Earth's interior involves searching for the model parameters which can match observed seismic data. We therefore require a forward model, or a way to compute the observations that we would expect to see if the earth had any given structure. Ideally, this forward model needs to encapsulate the full physics of wave propagation. If it does not, the seismograms predicted for each candidate structure will be inaccurate, and the results of the imaging process may be biased by this systematic error (for a full discussion, see Valentine & Trampert 2016). 2.1 The forward problem ‘Normal mode seismology’ represents one framework for constructing a seismological forward model, particularly suited to low-frequency, global-scale simulations. The theory incorporates all important physical effects, and can be used to produce seismograms that faithfully reproduce physical reality. However, the computational resources required are significant, rendering full calculations intractable for early studies. As a result, a variety of simplifications were developed. These can greatly reduce the computational costs, at the expense of a degradation in accuracy. A brief, self-contained account of the underlying theory is given in Appendix, emphasizing how and why these various approximations can be introduced. Here, we sketch some key points. 2.1.1 Reference eigenfunctions For any spherically symmetric reference model (e.g. PREM, Dziewonski & Anderson 1981), it is relatively straightforward to compute eigenfunctions $$\boldsymbol {\mathcal {S}}_k$$ and their associated eigenfrequencies ωk; each pair $$(\boldsymbol {\mathcal {S}}_k,\omega _k)$$ represents one ‘normal mode’ (e.g. Woodhouse 1988). Conventionally, four indices are used to describe the modes: k → (q, l, m, n), with q denoting the ‘class’ (principally ‘spheroidal’, q = S, or ‘toroidal’, q = T ); l, the angular degree, and m, the azimuthal order, identifying the spherical harmonic Ylm(θ, ϕ) which describes the lateral form of $$\boldsymbol {\mathcal {S}}_k$$; and n, the overtone number, identifying its radial behaviour. The eigenfrequency can be shown to depend only upon q, l and n. For each (q, l, n) combination, there are therefore 2l + 1 eigenfunctions that have a different spatial form, but share a common eigenfrequency. This degenerate set is often referred to as a ‘multiplet’, and identified by a label in the form nql(e.g. 0S2). The constituent eigenfunctions of each multiplet are then referred to as ‘singlets’. 2.1.2 Computing synthetic seismograms A central theme of normal mode seismology is that these eigenfunctions, obtained for a spherically symmetric reference model, may be used as a convenient set of basis functions for representing vector fields within the (laterally heterogeneous) Earth. In particular, the seismic waveforms s(x, t) can be expressed in the form (cf. eq. A5)   \begin{equation} \mathbf {s}(\mathbf {x},t) = \sum _k u_k(t)\boldsymbol {\mathcal {S}}_k(\mathbf {x}) \end{equation} (1)where $$\boldsymbol {\mathcal {S}}_k$$ is a vector field representing an eigenfunction of the reference model. There are an infinite number of eigenfunctions $$\boldsymbol {\mathcal {S}}_k$$, and so in principle this representation of the seismic wavefield requires an infinite summation. However, if seismograms are to be filtered so that they contain no component above a given frequency, ωmax, it is justifiable to neglect modes with eigenfrequency ωk ≫ ωmax. To make this concrete, we define a cut-off frequency, ωc ≥ ωmax, and work with the finite set of modes having ωk ≤ ωc. Then, the Fourier-transformed expansion coefficients $$\tilde{u}_k(\omega )$$ may be found by solving a system of linear equations (cf. eq. A9)   \begin{equation} \mathbf {M}(\omega ) \mathbf {\tilde{u}} = \mathbf {\tilde{q}} \end{equation} (2)Here, M(ω) is a frequency-dependent ‘coupling matrix’ which can be computed for any earth model presenting 3-D variations in density, elastic constants and attenuation, while the vector $$\mathbf {\tilde{q}}$$ is a representation of the seismic source (typically an earthquake). 2.1.3 Full-coupling The term ‘full-coupling’ implies that we solve this linear system (i.e. eq. 2) directly. If we do this, we obtain seismograms that are essentially physically complete. They incorporate only one fundamental approximation: truncation of the infinite series at frequency ωc. The first goal of this paper is to quantify the effect of this truncation upon synthetic seismic spectra, and to assess where ωc should lie in order for the truncation to have no appreciable impact upon calculations for frequencies below ωmax. 2.1.4 The self- and group-coupling approximations As ωc increases, the dimension of the matrix M increases, and thus the cost of solving eq. (2) grows rapidly. The elements of M are found by computing integrals of the form (cf. eq. A10)   \begin{equation} M_{ij}(\omega ) = \int _V \boldsymbol {\mathcal {S}}_i^*(\mathbf {x})\cdot \mathcal {A}(\omega )\boldsymbol {\mathcal {S}}_j(\mathbf {x})\,\mathrm{d}^3\mathbf {x} \end{equation} (3a)where $$\mathcal {A}(\omega )$$ is an integro-differential operator that depends upon frequency, and also upon the earth model. The elements tend to be relatively small unless $$\boldsymbol {\mathcal {S}}_i$$ and $$\boldsymbol {\mathcal {S}}_j$$ have similar eigenfrequencies. This leads to structure within M, linked to the degeneracy within each multiplet; it is dominated by the elements lying close to the diagonal. Motivated by this, the ‘self-coupling approximation’ modifies the definition of M, introducing an over-riding assertion that   \begin{equation} M_{ij}(\omega ) = 0 \,\,\textrm{if}\,\,\boldsymbol {\mathcal {S}}_i\ \textrm{and}\ \boldsymbol {\mathcal {S}}_j\,\,{\rm come\ from\ different\ multiplets.} \end{equation} (3b)This results in M taking a block-diagonal form, with each multiplet nql contributing a block of dimension (2l + 1) × (2l + 1). Solution of eq. (2) is then much easier, as it decouples into separate calculations for each multiplet. Of course, the resulting solution is only an approximation to the one that would be obtained if the full-coupling matrix were used. The ‘group-coupling’ approximation allows for limited interaction between multiplets that are adjacent in frequency. The spectrum is divided into ‘groups’, with each group containing one or more multiplets. Studies may differ in the groupings used. Then, eq. (3b) is relaxed slightly, so that instead   \begin{equation} M_{ij}(\omega )= 0 \,\,\textrm{unless}\,\,\boldsymbol {\mathcal {S}}_i\ {\rm and}\ \boldsymbol {\mathcal {S}}_j\,{\rm belong\, to\, the\, same\, group.} \end{equation} (3c) Again, this simplifies the solution of eq. (2), allowing separate computations for each group. Studies have shown that the resulting seismograms are more accurate than those of self-coupling, but errors remain significant (e.g. Deuss & Woodhouse 2001; Al-Attar et al.2012; Yang & Tromp 2015). 2.1.5 Normal modes of an aspherical earth model So far, we have outlined how to obtain seismograms in heterogeneous earth models. It is also possible to compute the normal modes of a general (non-spherically symmetric) earth model (see Appendix A4) and again, group- or self-coupling approximations may be adopted if required. We find that the modes are no longer degenerate: each singlet within a multiplet has a distinct eigenfrequency. This property is sometimes referred to as ‘splitting’ of the multiplet, relative to the spherically symmetric reference. Given that realistic earth models contain only weak lateral heterogeneity, the typical frequency shifts involved are small, and due to finite-length time series, we do not generally expect to be able to resolve isolated singlets in spectral observations. Nevertheless, aspherical structure manifests itself in the overall shape of each ‘peak’ in the spectrum of a seismogram. 2.2 From spectra to structure: the splitting functions Given the ability to compute synthetic spectra, one may contemplate formulating an inverse problem that aims to infer earth structure from observations. This could be implemented in a variety of ways, each with pros and cons. One popular approach involves the use of ‘splitting functions’. These were introduced within Woodhouse & Giardini (1985), and developed in detail by Giardini et al. (1987, 1988). These original papers are developed under the self-coupling assumption, where each multiplet is regarded as entirely isolated from every other. This was recognized to be a simplification of the underlying physics, and therefore an approximation, but it was necessary if calculations were to be feasible using the computational resources then available. Making this approximation, the contribution of a given multiplet κ = nql within seismograms is entirely characterized by the (2l + 1) × (2l + 1) diagonal sub-block in the coupling matrix arising from eqs. (3a & 3b), which we denote M(κκ). Similarly, the multiplet can be directly identified with a specific portion of an observed seismic spectrum. Thus, given a suite of observations of the relevant portion at different locations and from various seismic events, one can formulate an inverse problem to estimate the elements of M(κκ). Splitting functions provide a mechanism for parametrizing this inversion, designed to provide a straightforward link between M(κκ) and the underlying earth model. 2.2.1 Splitting functions for isolated groups of multiplets The splitting function framework was extended by Resovsky & Ritzwoller (1995) and others, allowing for cross-coupling between specific pairs (or small groups) of multiplets with similar frequencies as in group-coupling. For simplicity, we describe here the case where two multiplets, κ = nql and $$\kappa ^\prime = {}_{n^\prime }q^\prime _{l^\prime }$$, are allowed to interact with one another, but are assumed to be isolated from the remainder of the spectrum. Extension to larger groups, or restriction to the self-coupling case, follows straightforwardly. The isolated-group assumption implies that the only relevant part of M(ω) is the symmetric diagonal block   \begin{equation} \left(\begin{array}{@{}cc@{}}\mathbf {{M}}^{(\kappa \kappa )} & \mathbf {{M}}^{(\kappa \kappa ^\prime )}\\ {\mathbf {{M}}^{(\kappa ^\prime \kappa )}}&\mathbf {{M}}^{(\kappa ^\prime \kappa ^\prime )}\end{array}\right). \end{equation} (4)As discussed in Appendix A4, some additional assumptions allow us to introduce the splitting matrix, H, such that $$\mathbf {{M}}^{(\kappa \kappa ^\prime )} \approx \mathbf {H}^{(\kappa \kappa ^\prime )}-\omega ^2\mathbf {I}$$. We then parametrize this by introducing a set of complex-valued ‘splitting coefficients’ cst such that   \begin{equation} \mathbf {H}^{(\kappa \kappa ^\prime )}_{mm^\prime } = \omega _0^2\delta _{\kappa \kappa ^\prime } + \omega _0 \mathbf {W}_{mm^\prime }^{(\kappa \kappa ^\prime )}+\sum _{s=l - l^\prime }^{l+l^\prime }\sum _{t=-s}^s\gamma _{ll^\prime s}^{mm^\prime t} c_{st}^{(\kappa \kappa ^\prime )} \end{equation} (5)where the γ are defined in terms of spherical harmonic triple products,   \begin{equation} \gamma _{ll^\prime s}^{mm^\prime t} = \int _0^{2\pi }\int _0^\pi Y^*_{lm}(\theta ,\phi )Y_{l^\prime m^\prime }(\theta ,\phi )Y_{st}(\theta ,\phi )\sin \theta \,\mathrm{d}\theta \mathrm{d}\phi . \end{equation} (6)This integral can be expressed straightforwardly in terms of the Wigner-3j symbols, and vanishes for many combinations of spherical harmonics. Together, the coefficients describe a ‘splitting function’   \begin{equation} \eta ^{(\kappa \kappa ^\prime )}(\theta ,\phi ) = \sum _{s=l-l^\prime }^{l+l^\prime }\sum _{t=-s}^s c_{st}^{(\kappa \kappa ^\prime )} Y_{st}(\theta ,\phi ) \end{equation} (7)which may be regarded as a map depicting spatial variations in spectral properties. Given the symmetry of M, three splitting functions are required to fully characterize our isolated pair of multiplets: η(κκ), $$\eta ^{(\kappa ^\prime \kappa ^\prime )}$$ and $$\eta ^{(\kappa \kappa ^\prime )}$$. Following Woodhouse & Girnius (1982), the contribution of the two multiplets within each seismogram can be computed via diagonalization of eq. (4) (see eq. A13 and the associated discussion). Thus, it is possible to pose the inverse problem whereby the cst are determined from observed data; typically, the fiducial frequency ω0 and parameters describing attenuation (which we will not discuss further) are also refined within this process. The inversion is found to be weakly non-linear, and is typically tackled using a Bayesian formulation of the least-squares algorithm (as per Tarantola & Valette 1982). Invariably, the data is insufficient to uniquely constrain a solution, and it is therefore necessary to introduce prior information (i.e. regularization). Discussion of the extent to which this influences results may be found in numerous sources (e.g. Resovsky & Ritzwoller 1999; Kuo & Romanowicz 2002; Pachhai et al.2016). In the isolated-group approximation, the splitting coefficients can be straightforwardly related to structural parameters. For the purposes of illustration, we assume that the aspherical earth model specifies P-wave speed (α), S-wave speed (β), and density (ρ) perturbations at every point with respect to a spherically symmetric reference model. Furthermore, we assume that the model is laterally parametrized in terms of spherical harmonics to maximum degree L. Thus, the P-wave speed field may be expressed   \begin{equation} \alpha (r,\theta ,\phi ) = \sum _{l=0}^L \alpha _{lm}(r) Y_{lm}(\theta ,\phi ) \end{equation} (8)with identical expressions for β and ρ. By deriving the exact expression for $$\mathbf {H}^{(\kappa \kappa ^\prime )}_{mm^\prime }$$ (e.g. Woodhouse 1980), and then comparing this with eq. (5), it is possible to identify formulae for the kernels K such that   \begin{eqnarray} c^{(\kappa \kappa ^\prime )}_{st} &=& \int _0^a \alpha _{st}(r)K^{(\alpha ;\kappa \kappa ^\prime )}_{st}(r)+\beta _{st}(r) K_{st}^{(\beta ;\kappa \kappa ^\prime )}(r)\nonumber\\ &&+\,\, \rho _{st}(r)K^{(\rho ;\kappa \kappa ^\prime )}_{st}(r)\,\mathrm{d}r. \end{eqnarray} (9)If the aspherical model is specified in terms of a different set of parameters, equivalent expressions can be found. Thus, an earth model can be obtained by performing a linear inversion upon measured splitting coefficients. A particular attraction of this procedure is the possibility of selecting multiplets that are presumed to have sensitivity at particular depths, and hence to produce models focused on a given region of the Earth (e.g. Koelemeijer et al.2013). 2.2.2 Splitting functions in the real Earth All the foregoing discussion is self-consistent, subject to the fundamental assumption that the (groups of) multiplets are spectrally isolated. However, at least in the context of forward-modelling, it has been conclusively shown that this approach is insufficiently accurate for modern seismology, and instead high-quality seismic calculations must account for interactions throughout the spectrum. Various studies have shown that allowing for cross-coupling in splitting function analysis improves results (e.g. Resovsky & Ritzwoller 1995, 1998; Irving et al.2008; Deuss et al.2013). However, it is not possible to extend the splitting function approach to ever-wider coupling groups: the linearization described in Appendix A4 requires the coupling band-width to be small. If, instead, we choose to work with the full frequency-dependent coupling matrix M(ω), it is unclear whether a splitting-function–style parametrization can be usefully introduced. Where does this leave splitting-function studies? In effect, splitting functions should be regarded as a derived data type or ‘secondary observables’ (e.g. Cara & Lévêque 1987), intended to convey spectral information in a more manageable form. The measurement procedure should therefore be seen as a data-processing operation, or transformation applied to raw seismic data. The particular form of this transformation may be motivated by an out-dated assumption, but it remains a valid transformation to adopt. As such, published measurements of splitting functions may be taken at face value: they need not be regarded as somehow ‘approximate’. The splitting function is defined to be the result of a specified measurement procedure. However, if we are to adopt the isolated-multiplet assumption in the measurement procedure, eq. (9) provides only an approximate relationship between the measured splitting functions and the underlying earth structure. Alternatively, one may take the view that eq. (9) is the defining property of the splitting functions: in this case, the measurement procedure is deficient, since it does not properly account for physical effects in wave propagation. Regardless of the viewpoint adopted, ‘measured’ and ‘predicted’ splitting functions differ, due to approximations within the measurement procedure. Thus, the second key goal of this paper is to assess the impact of this approximation: can models produced using eq. (9) nevertheless be interpreted usefully, or are they too severely contaminated by systematic errors? If fully accurate earth models are to be produced from splitting function data, it is necessary to replace eq. (9) and properly account for multiplet interactions throughout the spectrum. Unfortunately, it is not clear that this can be achieved by simple modifications of the modal depth functions K(r). It seems likely that the most fruitful course would involve obtaining sensitivity kernels via adjoint techniques within a full-coupling framework, allowing the full range of physical effects to be accounted for. However, we do not pursue this idea further within the current paper. Indeed, given the capacity to compute exact sensitivity kernels for spectra using adjoint methods, one may question whether there is value in continuing to adopt the two-stage inversion procedure inherent to splitting function studies. 3 TRUNCATION ERRORS IN FULL-COUPLING SEISMOGRAMS First, we consider the problem of computing seismograms using full-coupling. As we have described, a theoretically complete expression for normal mode seismograms requires summation of an infinite set of seismograms. To make the computations tractable, we must truncate this modal expansion (eq. 1), and only include those modes with eigenfrequencies below some cut-off, ωc. This truncation will inevitably introduce errors into the computed seismograms. How do we choose ωc to ensure that synthetic spectra are accurate below some frequency ωmax? 3.1 The scale of acceptable error In order to address this question, we must somehow define the magnitude of error that we consider acceptable. One way to go about this is to consider the noise in the observed seismograms that we ultimately wish to fit. To allow meaningful analysis, the truncation error should be well below that noise. However, noise levels may vary considerably between different locations and experiments, limiting the generality of this definition. In addition, the earthquake source itself is also not perfectly well known, although it is apparent from inspection of eq. (2) that this uncertainty can be absorbed into the observational data uncertainty, by adding the corresponding variances. Alternatively, given that we typically wish to fit seismograms to infer Earth structure, another definition of an ‘acceptable’ truncation error is that the error should be much smaller than the signal we wish to infer. This is perhaps more straightforward, since it does not depend upon the vagaries of seismic noise. Of course, the signal itself must be above the noise level if we are to make meaningful inferences. We discuss both cases below and show that they can lead to different conclusions depending on the earth model. To quantify these ideas, we will only consider synthetic data calculated for various earth models. Of the three major seismic parameters (S-wave speed, P-wave speed & density), density is known to have the least influence on seismograms, and is a parameter of considerable importance for understanding global geodynamics. We therefore focus on the density signal when assessing whether truncation errors are ‘negligible’. 3.2 Observations using S20RTS We solved eq. (2) exactly, taking all 3-D effects of varying elastic constants and density into account, using the method described in Al-Attar et al. (2012), and computed the same set of seismograms several times varying ωc. In all cases, we used the Global CMT Catalog source mechanism for the 1994 Bolivian event (event code: 060994A) and a global distribution of 129 stations based on those within the Global Seismic Network and IRIS/IDA Seismic Network (see Fig. 1 ). As basis functions, we used eigenfunctions calculated in in the spherically symmetric model PREM (Dziewonski & Anderson 1981). The coupling matrix, explicitly defined in eq. (A10a), is constructed for the shear wave speed model S20RTS (Ritsema et al.1999), which provides values for δ ln Vs with respect to PREM, and where we added scaled compressional wave speed perturbations using δ ln Vp = 0.5 δ ln Vs and scaled density perturbations using δ ln ρ = 0.3 δ ln Vs. These scaling relations are commonly used and are obtained from mineral physics considerations (e.g. Karato 1993). On top of the mantle model, we implemented the simple crustal model from Woodhouse & Dziewonski (1984) and included rotation and ellipticity. In the following, we will consider vertical component spectra. Figure 1. View largeDownload slide Position of the stations (red stars) and the Bolivian event (beach ball) used in this study. Station GRFO is highlighted and shown by a blue star. Figure 1. View largeDownload slide Position of the stations (red stars) and the Bolivian event (beach ball) used in this study. Station GRFO is highlighted and shown by a blue star. We chose ωmax = 3 mHz and used ωc = 3,  4,  5,  and 6 mHz, increasing the number of coupled singlets (toroidals and spheroidals combined) from approximately 2000 to 13 000. We treated the results from ωc = 6 mHz as a set of ‘reference’ synthetic spectra, $$\mathbf {s}_\mathrm{ref}$$, and looked at amplitude spectra of the differences $$\mathbf {s}_{(\omega _c=3\,\mathrm{mHz})}-\mathbf {s}_\mathrm{ref}$$, etc., to provide an indication of the convergence of the full-coupling calculations (Fig. 2). These difference spectra represent the truncation error, and are referred to as such in the figure. It is clear that a truncation at ωc = ωmax can still have substantial errors, but they decrease rapidly as ωc is increased. As far as we know, most published studies that adopt a ‘full-coupling’ strategy have used ωc ≈ ωmax. Figure 2. View largeDownload slide Amplitude spectrum of vertical acceleration at station GRFO (black line) calculated by coupling all modes having frequencies less than ωc = 6 mHz. We refer to this as the ‘reference’ spectrum, sref. Also shown are amplitude spectra of the difference between this reference spectrum and those produced using lower values of ωc (blue lines), which we describe as the ‘truncation error’ for that cut-off. The red line shows the effect of removing any contribution from lateral density variations within the reference spectrum. Note that this is of broadly similar scale to the truncation error at ωc = 3 mHz. The green line shows the effect of perturbing the seismic source with three standard deviations on the reported source uncertainties. Figure 2. View largeDownload slide Amplitude spectrum of vertical acceleration at station GRFO (black line) calculated by coupling all modes having frequencies less than ωc = 6 mHz. We refer to this as the ‘reference’ spectrum, sref. Also shown are amplitude spectra of the difference between this reference spectrum and those produced using lower values of ωc (blue lines), which we describe as the ‘truncation error’ for that cut-off. The red line shows the effect of removing any contribution from lateral density variations within the reference spectrum. Note that this is of broadly similar scale to the truncation error at ωc = 3 mHz. The green line shows the effect of perturbing the seismic source with three standard deviations on the reported source uncertainties. To provide a sense of scale for this truncation error, we repeated the calculation using ωc = 6 mHz, but using a version of S20RTS where all lateral density heterogeneity has been omitted (i.e. with density structure as given by PREM). The amplitude spectrum of $$\mathbf {s}_{\bar{\rho }}-\mathbf {s}_\mathrm{ref}$$ is also shown in Fig. 2 and gives some indication of the magnitude of seismic signal that might be attributable to density structure. Perhaps surprisingly, we observe that this is similar to the magnitude of errors arising from truncation at ωc = ωmax. To estimate the effect of source uncertainties, we make use of those reported in the Global CMT Catalogue, although we note that these are likely to be an under-estimate (Valentine & Trampert 2012). We therefore implement a large source variation by perturbing each parameter by three times the corresponding reported standard deviation, with positive or negative perturbations chosen randomly, and again computed spectra. Given the other effects seen in Fig. 2, the source effect is very modest. We will therefore neglect its contribution in the discussion below. Rather than plotting similar figures for every station around the globe, we sought a method to assess the truncation errors more systematically. We defined a measure of the relative spectral misfit within any given frequency range (ω1, ω2)   \begin{equation} \xi (\mathbf {s},\omega _1,\omega _2) = \frac{\int _{\omega _1}^{\omega _2}\left|\mathbf {s}(\omega )-\mathbf {s}_\mathrm{ref}(\omega )\right|^2\,\mathrm{d}\omega }{\int _{\omega _1}^{\omega _2} \left|\mathbf {s}_\mathrm{ref}(\omega )\right|^2\,\mathrm{d}\omega }. \end{equation} (10)We then plot histograms of ξ evaluated for every station in our data sets (Fig. 3). Again, we compare the reference spectra to those truncated at ωc = 3, 4, and 5 mHz, and to those obtained by using ωc = 6 mHz, but omitting lateral density. If we wish to ensure that the truncation error is consistently smaller than the density signal, we require the corresponding histograms to not overlap. We clearly see that a truncation at 3 mHz will generally lead to errors of a similar scale to the signal anticipated from density. Even for a truncation at 4 mHz some overlap of the distributions occurs for frequencies close to ωmax. It therefore seems necessary to ensure that the truncation frequency is significantly above the maximum frequency of interest (i.e. ωc ≫ ωmax). These results are broadly in accordance with what we might expect based on the knowledge of the coupling-matrix elements. Generally, modes close together in frequency interact most strongly, so truncation errors ought to grow as we approach the cut-off frequency. Figure 3. View largeDownload slide Histograms showing relative spectral misfit, ξ (eq. 10), across all stations. Again, we compare reference spectra (ωc = 6 mHz) to those truncated at lower frequencies (ωc = 3, 4, 5 mHz), with all calculations being performed in the model S20RTS (where we incorporated Vp and density heterogeneity as described in the main text). Results from comparing the same reference spectra to those calculated in a version of S20RTS without any lateral density heterogeneity are also shown. Figure 3. View largeDownload slide Histograms showing relative spectral misfit, ξ (eq. 10), across all stations. Again, we compare reference spectra (ωc = 6 mHz) to those truncated at lower frequencies (ωc = 3, 4, 5 mHz), with all calculations being performed in the model S20RTS (where we incorporated Vp and density heterogeneity as described in the main text). Results from comparing the same reference spectra to those calculated in a version of S20RTS without any lateral density heterogeneity are also shown. To quantify the truncation effects further, we defined a cumulative misfit for the entire synthetic data set   \begin{equation} \psi (\omega ) = \frac{\sum _ i \int _0^\omega |\mathbf {s}^{(i)}(\omega ^\prime )-\mathbf {s}_\mathrm{ref}^{(i)}(\omega ^\prime )|^2\,\mathrm{d}\omega ^\prime }{\sum _i \int _0^\omega |\mathbf {s}_\mathrm{ref}^{(i)}(\omega ^\prime )|^2\,\mathrm{d}\omega ^\prime } \end{equation} (11)where the summation is over all the stations in the data set. We then plot this quantity as a function of frequency for the various cut-off frequencies (Fig. 4), and for the data set where density heterogeneity has been omitted. There are several ways of reading this figure. First, focusing on the curves for the various choices of ωc: suppose we want a precision of (say) 0.1 per cent of the total energy in the signal (this choice should be guided by knowledge of the uncertainties in observed seismic spectra and the seismic source). In this case, truncation at ωc = 3 mHz will produces seismograms accurate up to ωmax = 2 mHz (as shown by a horizontal line on the figure). For a truncation at ωc = 4 mHz, seismograms are accurate up to 3 mHz within the same precision. It appears that we roughly need coupling of 1 mHz higher than the highest frequency of interest in the seismogram. Alternatively, if we are interested in imaging Earth's density structure, we should require truncation errors to be at least an order of magnitude smaller than the density signal in the spectra, and we see roughly that the cut-off frequency should be at 4 mHz for any frequency ωmax ≤ 3 mHz. These values are based on an assumption that Earth's heterogeneity is similar to that present in our implementation of S20RTS. The reader may choose her/his own desired precision and adapt the conclusions accordingly, by direct reading of Fig. 4. Figure 4. View largeDownload slide Cumulative misfits defined by eq. (11), relative to reference spectra computed in S20RTS, for three different cut-off frequencies and for the omission of lateral density heterogeneity. The horizontal grey line denotes a misfit of 0.1  per  cent, as discussed in the main text. Figure 4. View largeDownload slide Cumulative misfits defined by eq. (11), relative to reference spectra computed in S20RTS, for three different cut-off frequencies and for the omission of lateral density heterogeneity. The horizontal grey line denotes a misfit of 0.1  per  cent, as discussed in the main text. 3.3 Model-dependence of the truncation error We expect results to vary depending on the model, and in particular on the power spectrum of heterogeneity. Coupling effects are somewhat analogous to scattering, and it is well-known that short-wavelength heterogeneity promotes multiple-scattering effects. We might reasonably expect to observe similar effects in coupling calculations, with stronger short-wavelength structure necessitating relatively higher values of ωc. S20RTS has modest amplitudes and a ‘red’ spectrum, dominated by long-wavelength structure. The model resolution is far below what the parametrization of the model would allow, ranging from 2000–4000 km horizontally to 250–750 km vertically (Ritsema et al.2004). It is therefore important to investigate how this is dictating the findings in the previous subsection. To provide some insight into this question, we generated a random model using the same parametrization as S20RTS—a spherical harmonic expansion to degree 20 laterally and 21 splines vertically. For δ ln Vs we drew random numbers uniformly in the range [−0.01, 0.01] giving an rms-amplitude of about 12 per cent at all depths. We continue to assume that Vp and ρ are scaled versions of Vs, as with S20RTS. The rms-amplitude of this random model is between 5 and 38 times stronger than S20RTS depending on depth. This produces a model of true degree-20 lateral resolution, a white spectrum, and much higher vertical resolution. We repeated all the full-coupling calculations already described, using this model. Again, we produce histograms of relative spectral misfit (Fig. 5) and cumulative misfit (Fig. 6 ). From Fig. 5, we see that now the truncation effects are much stronger and overall we need to couple more modes than was necessary in S20RTS. In other words, more power in the structural model, and especially at short-wavelengths, necessitates higher values of ωc. It is also obvious that the contribution of density to the seismograms is much lower in the random model than in S20RTS, despite the fact that the rms-amplitude of the model is much larger. This is because density imprints itself differently into the seismograms, compared to the model with a red spectrum. Figure 5. View largeDownload slide Histograms of relative spectral misfit across all stations defined by eq. (10) for a random model. Figure 5. View largeDownload slide Histograms of relative spectral misfit across all stations defined by eq. (10) for a random model. Figure 6. View largeDownload slide Cumulative misfits defined by eq. (11) for a randomly generated, spectrally white earth model. Compare Fig. 4. Figure 6. View largeDownload slide Cumulative misfits defined by eq. (11) for a randomly generated, spectrally white earth model. Compare Fig. 4. Clearly the relation between ωc, ωmax and the earth model is not straightforward. Looking at cumulative misfits (Fig. 6), we see that white models appear to yield larger truncation errors than their red counterparts, as expected. Again, considering the data sets that include density heterogeneity, we can imagine seeking a precision of 0.1 per cent of the total energy in the signal. A truncation at ωc = 3 mHz now only produces seismograms accurate up to 1 mHz. For a truncation at ωc = 4 mHz, seismograms are accurate up to around 2 mHz within the same precision. It appears that we now need coupling to be at least 2 mHz higher than the highest frequency of interest in the seismograms. If, however, we define the acceptable truncation error as being an order of magnitude smaller than the signal arising from density heterogeneity, we now see that coupling up to 5 mHz is required for ωmax ≤ 3 mHz. Again it is straightforward to adapt the conclusions for any other desired precision. 3.4 Summary The key conclusions we draw from these experiments may be summarized: Accurate full-coupling spectra require the cut-off frequency ωc to be significantly above the maximum frequency present in the seismograms, ωmax. To our knowledge, this criterion has not been met in previously published studies based on full-coupling. Accurate seismograms up to 3 mHz in earth-like models (having a red spectrum) require truncation at no less than 4 mHz. Generally truncation errors are stronger at higher frequencies, as might be expected. Effects seem more pronounced as models gain more power, especially at shorter wavelengths. In models with white spectra, accurate seismograms up to 3 mHz require a truncation level of at least 5 mHz. Overall, these observations are all broadly consistent with expectations based on theoretical arguments. However, the effects are now quantified, and we have obtained rules-of-thumb that can be applied to ensure acceptable accuracy when performing calculations. 4 THE INTERPRETABILITY OF SPLITTING FUNCTION MEASUREMENTS We now turn to consider splitting functions, and the extent to which they can provide useful constraints upon Earth structure. Catalogues of splitting function measurements (e.g. Giardini et al.1987; He & Tromp 1996; Resovsky & Ritzwoller 1998; Durek & Romanowicz. 1999; Masters et al.2000; Deuss et al.2011, 2013) still represent the main ‘derived’ data set for long-period seismology. Measured centre frequencies (i.e. ω0 in eq. 5) are taken as representative of the Earth's spherical average structure (e.g. de Wit et al.2014) and splitting function coefficients feature as a constraint in many current global tomography models (e.g. Ritsema et al.2011; Moulik & Ekström 2016). Splitting functions have also been used for targeted studies investigating specific aspects of Earth structure, since one can examine the kernels from eq. (9) and identify modes which have the desired sensitivity (see Fig. 7). This has particularly been exploited to investigate deep-earth structure, including inner core studies (beginning with Woodhouse et al.1986), targeted lower mantle studies (e.g. Koelemeijer et al.2013, 2017) and density inferences (e.g. Ishii & Tromp 1999; Resovsky & Trampert 2003). Figure 7. View largeDownload slide Self-coupling sensitivity kernels $$K_{00}^{(\kappa \kappa )}$$ (as in eq. 9) corresponding to the angular order s = 0 splitting coefficient for our selected modes. The solid red line denotes Vp sensitivity, the solid black line Vs sensitivity and the dotted black line gives density sensitivity. The modes are ordered somewhat arbitrarily, with a general trend from mostly upper-mantle–sensitive to mostly lower-mantle–sensitive. Figure 7. View largeDownload slide Self-coupling sensitivity kernels $$K_{00}^{(\kappa \kappa )}$$ (as in eq. 9) corresponding to the angular order s = 0 splitting coefficient for our selected modes. The solid red line denotes Vp sensitivity, the solid black line Vs sensitivity and the dotted black line gives density sensitivity. The modes are ordered somewhat arbitrarily, with a general trend from mostly upper-mantle–sensitive to mostly lower-mantle–sensitive. It is common for splitting function catalogues to publish uncertainty estimates on the recovered splitting coefficients and centre-frequencies. In principle, these can be obtained directly from the Bayesian inversion framework (e.g. Tarantola & Valette 1982). However, it is not clear how to define a meaningful prior for the splitting function inversion, and in any case such analysis would neglect the inherent non-linearity of the measurement process. This point has previously been made by Resovsky & Ritzwoller (1998). After a careful analysis of all potential contributions to the uncertainty, they estimated combined uncertainties using regressions performed on synthetic data with simulated theoretical errors and data uncertainties. More recently, it has been common to use some form of boot-strapping procedure to investigate the degree of constraint upon results (e.g. Deuss et al.2013). This uncertainty estimate therefore reflects any apparent inconsistency between different parts of the data set: it may not reflect any systematic effects inherent to the measurement procedure or its implementation. Various studies have highlighted that considering interactions between adjacent modes can alter splitting-function results significantly (e.g. Resovsky & Ritzwoller 1995; Irving et al.2008) and group-coupling, where appropriate, is now the norm for various frequency bands (e.g. Deuss et al.2013). Studies have also investigated the extent to which splitting function inversions are well-constrained, typically focusing on aspects of the linearized inversion scheme, such as regularization, parametrization and starting model (e.g. Resovsky & Ritzwoller 1999; Kuo & Romanowicz 2002; Pachhai et al.2016). However, we are not aware of any previous studies that have taken a more holistic approach, to investigate whether splitting functions provide a meaningful representation of Earth structure that can be interpreted by eq. (9). As already discussed, a potential issue arises because the measurement procedure is inherently based upon isolated-group approximations (either self- or group-coupling), although their interpretation (via eq. 9) is not. Our results in Section 3, and those of earlier work (e.g. Deuss & Woodhouse 2001; Al-Attar et al.2012; Yang & Tromp 2015), all indicate that isolated-group seismograms provide a poor approximation to the true physics of wave propagation in the earth. It therefore seems distinctly possible that structural models estimated via isolated-group synthetics may inherit these inaccuracies. We address this by performing a straightforward synthetic test. For any earth model, we compute ‘predicted’ self-coupling splitting function coefficients using eq. (9). We also compute accurate synthetic seismograms, using full-coupling theory and the results of the previous section. From these seismograms, we then obtain ‘measured’ splitting functions $$c^{(\kappa \kappa )}_{st}$$, following the same procedure as is used for real data. We then assess the agreement between measured and predicted splitting function values. Ideally, they should agree to within the stated uncertainty for the measurement procedure. If this is the case, we can conclude that the measured splitting functions are interpretable using eq. (9), and the measurement uncertainties can be propagated through this expression to provide uncertainties on structural parameters. If measured and predicted splitting functions differ by more than the measurement uncertainty—but show no evidence of any systematic bias—it may be reasonable to continue to apply eq. (9) for interpretation using increased uncertainty levels to reflect the theory error (as per Tarantola & Valette (1982)). However, if systematic biases are found, no meaningful interpretation can be made. 4.1 Measurement of synthetic splitting functions As stated in the detailed analysis of Resovsky & Ritzwoller (1998), splitting function measurements depend on theoretical errors, regularization, noise in the data and event-station distributions. The final splitting function depends non-linearly on all of these. We want to quantify as precisely as possible the theoretical errors, and so we need to specify the other contributions as accurately as possible. To do this, all details matter. We therefore generated a synthetic full-coupling data set and measured splitting functions following the approach of Deuss et al. (2013). We considered these new spectra to be the observed data, referred to as above as $$\mathbf {s}_\mathrm{ref}$$. We will only focus on the self-coupling splitting function η(κκ) (eq. 7), although measurements in certain parts of the spectrum employed group-coupling. We used the same event and station distributions and method (forward theory, cut-off frequencies and regularization) as those in Deuss et al. (2013). The measured splitting coefficients $$c^{(\kappa \kappa )}_{st}$$ can then directly be compared to the expected coefficients for model S20RTS (as obtained from eq. 9). Two factors contribute to any difference between measured and expected coefficients. First, differences may be the result of the approximations in the measurement procedure: this is the effect we wish to address. However, they may also arise as a result of the regularization applied within the non-linear optimization procedure. To separate the two contributions, we also calculated S20RTS seismograms strictly using the self-coupling approximation, and again measured splitting functions from these. In this case, any difference between measured and observed coefficients can be attributed solely to regularization, as data and measurement procedure now contain the same physics. Assuming that all contributions are Gaussian distributed, the difference between the covariances corresponding to the two sets of splitting function then isolates the theoretical error. We focused on a selection of modes, which were also measured by Deuss et al. (2013), namely: 0S2, 0S5, 0S6, 0S7, 1S4, 1S7, 1S8, 1S9, 1S10, 2S6, 2S12, 2S13, 5S3. Fig. 7 shows some of the kernels K, which we used in eq. (9) to predict the splitting functions we ought to retrieve if the measurements were perfect. Some of the chosen modes have sensitivity mainly in the lower mantle, some mainly in the upper mantle, and some are sensitive throughout the mantle. Our discussion will concentrate on coefficients of angular order 2 and 4, which carry most of the splitting function energy for the modes we selected to analyse. This will also avoid the problem of neglected higher degree structure as the measurements are done to higher orders (Resovsky & Ritzwoller 1998). All our synthetic experiments are noise-free because we want to isolate the theoretical error without a particular noise model (which would be difficult to choose objectively) propagating non-linearly into the results via the regularization. Assuming that splitting function measurements are secondary data to be taken at face value, we will put our synthetic measurements into context by using actual uncertainties from Deuss et al. (2013) obtained by boot-strapping real earthquake data and thought to represent only the data noise present in the Deuss catalogue. 4.2 Observations from synthetic splitting functions The coefficients of the splitting function measurements for mode 0S6, a lower mantle mode, are depicted in Fig. 8. When self-coupling seismograms are used, the measured coefficients are indistinguishable from those calculated for S20RTS via the kernels, suggesting that the regularization used for this mode is not affecting the measurement. We therefore only show the calculated coefficients in the figure, but not those measured from self-coupling seismograms. There is however a difference with those measured using the full-coupling spectra as data. To calibrate this difference against the density signal, we also computed predicted splitting coefficients for a version of S20RTS without lateral density heterogeneity. While there are differences between the three sets of splitting coefficients, these are small compared to the magnitude of the coefficients. It is therefore instructive to plot differences. Figure 8. View largeDownload slide Observed splitting coefficients for mode 0S6 for s = 2 and s = 4 (top panels). The full black circles represent measurements made on synthetic spectra obtained using full-coupling. The red open circles are predicted splitting coefficients calculated using eq. (9), and correspond to those you would expect if the measurements were perfect. The grey diamonds represent splitting coefficients calculated for mantle model S20RTS without lateral density variations. Bottom panels: difference between those measured from full-coupling spectra and predicted coefficients (black circles) and difference between predicted coefficients without and with density (grey diamonds). The vertical bars correspond to error bars inferred by Deuss et al. (2013) for this mode. Figure 8. View largeDownload slide Observed splitting coefficients for mode 0S6 for s = 2 and s = 4 (top panels). The full black circles represent measurements made on synthetic spectra obtained using full-coupling. The red open circles are predicted splitting coefficients calculated using eq. (9), and correspond to those you would expect if the measurements were perfect. The grey diamonds represent splitting coefficients calculated for mantle model S20RTS without lateral density variations. Bottom panels: difference between those measured from full-coupling spectra and predicted coefficients (black circles) and difference between predicted coefficients without and with density (grey diamonds). The vertical bars correspond to error bars inferred by Deuss et al. (2013) for this mode. In the bottom panels of Fig. 8 we plot the difference between the splitting functions predicted for S20RTS, and those measured from full-coupling synthetics (black line). Each coefficient is plotted with error-bars, representing the uncertainty reported for this mode by Deuss et al. (2013). We see that most coefficients touch the zero-line when the error-bars are taken into account (certainly when 2 standard deviations or the 95 per cent confidence level is considered), implying that predicted and measured splitting functions agree to within the measurement uncertainty. Thus, for this mode, any systematic errors that arise from the use of self-coupling are negligible compared to the effect of noise within the data. We also plot the difference between splitting functions predicted in S20RTS, and our version without density (grey line). Again, we plot error-bars representing the uncertainty reported by Deuss et al. (2013), and find that many of these also touch the zero-line and certainly within 2 standard deviations. This suggests that the signal from density lies at, or below the noise level of measured splitting functions. As a result, it is not straightforward to extract meaningful information about density from splitting functions for this mode. In addition, we observe that the signal attributable to lateral density heterogeneity appears to be similar in magnitude to the difference between observed and predicted splitting coefficients. Thus, even if noise levels could be reduced significantly, it would appear that current splitting function measurement procedures are insufficiently accurate to make use of the density signal. Similar results are seen in Fig. 9, which shows results for the upper mantle mode 2S12. In this example, we see differences which are more significant relative to uncertainties. Errors due to the use of self-coupling in the measurement procedure are now markedly greater than the quoted uncertainty, which accounts only for effects due to noise in the data, even within 2 standard deviations for some coefficients. This is likely to affect any interpretation based on eq. 9. Although lateral density heterogeneity has a signal that is now also substantially greater than the stated noise level, it remains comparable in size to the theoretical error in splitting function observations. Again, it seems unlikely that meaningful density inferences can be made from this mode using standard methods. These two examples suggest that both density signal, and measurement error, become more significant for modes with shallow sensitivity. To see if this holds in general, we must examine more modes. Figure 9. View largeDownload slide As Fig. 8, but for mode 2S12. Figure 9. View largeDownload slide As Fig. 8, but for mode 2S12. Comparing the full suite of splitting function coefficients measured from full-coupling synthetics to those predicted using eq. (9), we find a very strong correlation (with an R-value of 0.99), with no obvious evidence of any systematic trends. Splitting functions depend non-linearly on the spectra, and a more sophisticated error analysis is required before the possibility of any systematic effects can be fully eliminated. However, our results suggest that a reasonable starting point may be to assume that theoretical measurement errors can be modelled by a zero-mean Gaussian distribution, allowing them to be straightforwardly incorporated into measurement uncertainties. To interpret our observations, it is then most instructive to show the reduced χ2 misfits for all measured coefficients of a given splitting function with respect to the coefficients predicted for S20RTS using eq. (9)   \begin{equation} \chi ^2 = \frac{1}{N} \sum _{s,t} \frac{\left(c^\mathrm{meas.}_{st}-c^\mathrm{pred.}_{st}\right)^2}{\sigma ^2_{st}} \end{equation} (12)where σst are the uncertainties from Deuss et al. (2013) for the given mode and N represents the number of coefficients in the corresponding splitting function. This statistical measure allows results to be summarized for easy interpretation; χ2 = 1 indicates that the difference between measured and predicted values is, on average, of similar size to the uncertainties. First, we consider the splitting function measurements made upon synthetic seismograms computed using full-coupling (see Fig. 10 ). For most modes, the value of χ2 is well above one, indicating that observations differ from predicted values by an amount significantly greater than currently published data uncertainties for splitting coefficients. In other words, the Deuss catalogue underestimates the true uncertainty on splitting coefficients. This in turn may impact upon the confidence that can be placed in models derived partly or fully from these splitting function data. Figure 10. View largeDownload slide Reduced χ2 misfit defined by eq. (12) for measurements on full- (black squares) and self-coupling (open squares) seismograms using all parameters as determined by Deuss et al. (2013). Figure 10. View largeDownload slide Reduced χ2 misfit defined by eq. (12) for measurements on full- (black squares) and self-coupling (open squares) seismograms using all parameters as determined by Deuss et al. (2013). As has already been discussed, discrepancies between measured and predicted splitting coefficients may arise due to reliance upon approximations in the measurement procedure, or as a consequence of regularization. To compare these effects, Fig. 10 also shows χ2 for splitting functions measured from synthetic seismograms generated using self-coupling. In this case, seismogram generation and splitting function measurement employ a consistent theoretical framework. Thus, any differences here arise solely from regularization. For several modes, we see that χ2 obtained from self-coupling synthetics is similar to that for full-coupling synthetics. This suggests that the discrepancy between prediction and observation in these cases is mostly attributable to regularization. It is possible that reducing the damping applied during the measurement procedure would improve the agreement between predictions and observations. However, this might also be expected to increase the sensitivity to data noise, resulting in increased error-bars on the recovered coefficients (Resovsky & Ritzwoller 1998). It is therefore not obvious that a net benefit would be gained; re-analysis of real data to investigate this rigorously is beyond the scope of the present study. Assuming that all contributions to the uncertainties are Gaussian distributed, covariances will simply add up, in other words the total posterior covariance in actually measured coefficients may be expressed as $$\mathbf {C}_\mathrm{tot}=\mathbf {C}_\mathrm{th}+\mathbf {C}_\mathrm{reg}+\mathbf {C}_\mathrm{d}$$, where $$\mathbf {C}_\mathrm{th}$$ is the theoretical uncertainty we wish to quantify, $$\mathbf {C}_\mathrm{reg}$$ is due to regularization in the regression and $$\mathbf {C}_\mathrm{d}=\sigma ^2_d \mathbf {I}$$ corresponds to data uncertainty including the influence of the event-station distribution. This decomposition is similar to that used in Resovsky & Ritzwoller (1998). The uncertainties quoted in Deuss et al. (2013) were determined using boot-strapping: randomly eliminating stations, and also whole events. This was done for a fixed regularization. Thus, these uncertainties therefore represent the data covariance, $$\mathbf {C}_\mathrm{d}$$, including the particular event-station distribution used to create the Deuss catalogue. Our measurements based on self-coupling seismograms (see Fig. 10) do not reflect the calculated splitting functions perfectly. Because these seismograms are noise free, the discrepancies may be a measure of the regularization noise, which can be expressed as a multiple of the data noise, $$\mathbf {C}_\mathrm{reg}=\alpha \sigma ^2_d \mathbf {I}$$. Values of α can be directly read from Fig. 10. In the same figure we also display the result for measurements on full-coupling seismograms. These represent the discrepancies due to regularization and the theoretical error combined, which correspond then to $$\mathbf {C}_\mathrm{th}+\mathbf {C}_\mathrm{reg}=\beta \sigma ^2_d \mathbf {I}$$. Hence β − α is a measure of $$\mathbf {C}_\mathrm{th}$$ expressed in units of the covariances in the Deuss catalogue. We see that $$\mathbf {C}_\mathrm{th}$$ and $$\mathbf {C}_\mathrm{reg}$$ changes for every mode, but as a quantitative indication: for the sample of modes we considered, on average β − α = 1.3 and α = 2.0. Thus for the Deuss catalogue, the theoretical measurement error, is on average 1.3 times the quoted variance and the regularization error is on average 2 times the quoted variance, resulting in $$\mathbf {C}_\mathrm{tot}=(1+1.3+2.0)\sigma ^2_d\mathbf {I}=4.3\sigma ^2_d\mathbf {I}$$. We therefore recommend users should multiply the published uncertainty in the Deuss catalogue by a factor of 2 to use a more representative number for the true measurement uncertainty. While the quoted numbers reflect the contribution to uncertainty in the Deuss catalogue, can they be used for other catalogues? The answer is no, because regularization errors and the data uncertainty translate differently into coefficient uncertainty, depending intimately on all the details on the measurement procedure (event-station distribution, regularization, number of iterations, cut-off frequencies, …). Because we do not have this information available for other catalogues, we cannot give detailed results for those. We can however say Resovsky & Ritzwoller (1998) made a real effort to quantify the true uncertainties given the computational restrictions of 30 years ago. To put these results into perspective, it is instructive to look at potential information on the Earth's structure contained in the splitting functions. We have already described how we generate a version of S20RTS omitting any lateral density heterogeneity, and use this to estimate the magnitude of signal attributable to density. Similarly, we can generate seismograms in versions of S20RTS that assume there is no lateral Vs heterogeneity (having only lateral variations in Vp and ρ), or no lateral Vp structure (only Vs and ρ). In each case, we can compute expected splitting coefficients using eq. (9), and obtain χ2 values as in eq. (12). When we plot these in Fig. 11, we see that the signal of Vs in the splitting functions is 2-3 orders of magnitude larger than the observational uncertainties for most modes, while that of Vp is on average one order of magnitude larger and ρ mostly of the same order as the uncertainties. This is interesting in itself, as it shows that most modes do not contain any significant information on density, regardless of whether we can accurately measure their splitting function or not. The reason for this might be that wave speeds mostly affect the phase of a seismic spectrum, whereas density mostly affects its amplitude. When we further consider our findings that the measured data variances should be multiplied by a factor of 4.3 to represent the total measurement errors, we find that on average the reduced χ2 will be around 100 for δ ln Vs, around 4 for δ ln Vp and around 1 for δ ln ρ. Figure 11. View largeDownload slide Reduced χ2 misfit defined by eq. (12) for model S20RTS, except with δ ln Vs = 0 (black squares), δ ln Vp = 0 (open squares) and δ ln ρ = 0 (grey diamonds). This provides an indication of the average strength of the respective fields within the splitting functions for the 13 modes analysed. Figure 11. View largeDownload slide Reduced χ2 misfit defined by eq. (12) for model S20RTS, except with δ ln Vs = 0 (black squares), δ ln Vp = 0 (open squares) and δ ln ρ = 0 (grey diamonds). This provides an indication of the average strength of the respective fields within the splitting functions for the 13 modes analysed. Most studies infer mantle structure using these splitting function in a least-squares inversion (e.g. Ritsema et al.1999; Ishii & Tromp 1999; Ritsema et al.2011; Moulik & Ekström 2016; Koelemeijer et al.2017). Given that the Vs signal is well above the uncertainty of the splitting functions, shear wave mantle models using such data should be fairly robust. However, because the density signal is of the level of the measurement uncertainty, its robustness was justifiably questioned in the past (e.g. Resovsky & Ritzwoller 1999; Kuo & Romanowicz 2002). Only by moving away from standard techniques and using full sampling, can density be inferred with a meaningful uncertainty (Resovsky & Trampert 2003). The latter study showed that the density model had on average an uncertainty of about 50 per cent of its magnitude, or 100 per cent within a 95 per cent confidence level. This generated lively debates, and for instance, the question of buoyant versus heavy large low-shear wave provinces is still waiting to be settled, although the heavy side (Ishii & Tromp 1999; Resovsky & Trampert 2003; Trampert et al.2004; Mosca et al.2012; Moulik & Ekström 2016) got a recent boost with an independent study inverting tidal data (Lau et al.2017). Given that the Deuss catalogue uses far more seismograms to infer splitting functions than any previous catalogue, and we now have a good representations of the true measurement uncertainty, constructing a mantle model using a sampling algorithm is well worth repeating. Alternatively, direct spectra inversions using full coupling theory, where density has a higher signal-to-noise ratio, are becoming computationally feasible. 4.3 Summary Again, our results can be summarized as: We confirm that there is a significant theoretical measurement error in currently published splitting functions. There is no obvious bias, but this theoretical error can, to a first approximation, be assumed to be random noise. There is also a significant contribution from regularization to the uncertainty in splitting functions. Taking the total uncertainty into account, currently published splitting functions contain information on Vs with a signal-to-noise ratio of about 100, information on Vp with a signal-to-noise ratio of about 4, and information on ρ with a signal-to-noise ratio of about 1. Since the density signal is small in currently published splitting functions, advanced techniques should be employed to infer density structure. The signal from Vs is very strong, and its inference is unlikely to be significantly affected by the increased measurement uncertainty in currently published splitting functions. Thus, we do not believe our results should affect published models of Vs structure that may be derived these splitting function measurements. These results are specific to the Deuss catalogue (Deuss et al.2013). Using other catalogues requires repeating the work in this study if authors wish to infer density, the exception possibly being the Colorado catalogue (Resovsky & Ritzwoller 1998). We also remark that the density signal may be more effectively isolated by working directly upon seismic spectra, where the observational uncertainty is expected to be smaller: the results of Section 3 suggest that this should be feasible. 5 CONCLUSIONS Previous studies have shown that self- and group-coupling approximations ought to be regarded as obsolete in the context of computing seismograms in 3-D Earth models. By systematically increasing the number of modes in full-coupling calculations, we have demonstrated that accurate full-coupling requires inclusion of modes significantly beyond the highest frequency range of interest in the seismograms. For seismograms to be accurate up to 3 mHz in typical models, coupling to 5 mHz is preferable; 4 mHz may be sufficient for some applications if the Earth has a structure similar to S20RTS. Our experimental evidence appears to corroborate theoretical expectations that truncation errors become larger as short-wavelength structure becomes more dominant. Self-coupling splitting functions, often used in the construction of 3-D Earth models, are themselves inherently based on the self- or group-coupling approximations, and may therefore contain systematic measurement errors. By measuring self-coupling splitting functions on exact full-coupling seismograms, we showed that the recovered splitting function does not match theoretical predictions. There is also a significant regularization contribution to the recovered splitting functions. There is no visible apparent bias from any of these contributions, but recently published (Deuss et al.2013) standard deviations of data uncertainties for splitting coefficients need to be increased by a factor of around 2 to represent the total error within the measurement process. In that case, our results show that the information on the 3-D density structure contained in the measured splitting functions has the same magnitude as the total uncertainty. The information on the variation of compressional-wave speed is about 4 times that of the uncertainty and the information on the shear-wave speed has a signal-to-noise ratio of about 100. To go beyond a statistical treatment of the theoretical error, one route forward may be to simply avoid the use of splitting functions: with modern computational resources, it may be feasible to infer structural information directly from spectra via full-coupling calculations. Alternatively, it may also be possible to continue to use splitting functions, measured using self- or group-coupling approximations, and instead replace the kernels within eq. (9) by those that account for these measuring approximations correctly (e.g. via adjoint theory). However, until such studies are undertaken, we must assume markedly larger uncertainties than reported for measured splitting coefficients. Given this, our modelling suggests that Earth's density signal lies at or below the presumed noise level. In consequence, there is no prospect of recovering robust estimates of density structure. It therefore follows that currently published density models ought to be interpreted with great caution. ACKNOWLEDGEMENTS We are grateful for constructive reviews by Philippe Lognonné and Joseph Resovsky. The research leading to these results has received funding from the European Research Council (ERC) under the European Union's Seventh Framework Programme (FP/2007-2013) grant agreement number 320639 (iGEO) and under the European Union's Horizon 2020 research and innovation programme grant agreement number 681535 (ATUNE). The code used to perform the full-coupling calculations makes extensive use of routines developed by John Woodhouse, and we thank him for allowing us to use them. REFERENCES Al-Attar D., Crawford O., 2016. Particle relabelling transformations in elastodynamics, Geophys. J. Int.  205 575– 593. https://doi.org/10.1093/gji/ggw032 Google Scholar CrossRef Search ADS   Al-Attar D., Woodhouse J., Deuss A., 2012. Calculation of normal mode spectra in laterally heterogeneous earth models using an iterative direct solution method, Geophys. J. Int.  189 1038– 1046. https://doi.org/10.1111/j.1365-246X.2012.05406.x Google Scholar CrossRef Search ADS   Cara M., Lévêque J., 1987. Waveform inversion using secondary observables, Geophys. Res. Lett.  14 1046– 1049. https://doi.org/10.1029/GL014i010p01046 Google Scholar CrossRef Search ADS   Dahlen F., 1968. The normal modes of a rotating, elliptical Earth, Geophys. J. R. astr. Soc.  16 329– 367. https://doi.org/10.1111/j.1365-246X.1968.tb00229.x Google Scholar CrossRef Search ADS   Dahlen F., 1969. The normal modes of a rotating, elliptical Earth—II. Near-resonance multiplet coupling, Geophys. J. R. astr. Soc.  18 397– 436. https://doi.org/10.1111/j.1365-246X.1969.tb03576.x Google Scholar CrossRef Search ADS   Dahlen T., Tromp J., 1998. Theoretical Global Seismology  Princeton University Press. de Wit R., Käufl P., Valentine A., Trampert J., 2014. Bayesian inversion of free oscillations for Earth's radial (an)elastic structure, Phys. Earth planet. Inter.  237 1– 17. https://doi.org/10.1016/j.pepi.2014.09.004 Google Scholar CrossRef Search ADS   Deuss A., Woodhouse J., 2001. Theoretical free-oscillation spectra: the importance of wide band coupling, Geophys. J. Int.  146 833– 842. https://doi.org/10.1046/j.1365-246X.2001.00502.x Google Scholar CrossRef Search ADS   Deuss A., Woodhouse J., 2004. Iteration method to determine the eigenvalues and eigenvectors of a target multiplet including full mode coupling, Geophys. J. Int.  159 326– 332. https://doi.org/10.1111/j.1365-246X.2004.02399.x Google Scholar CrossRef Search ADS   Deuss A., Ritsema J., van Heijst H., 2011. Splitting function measurements for Earth's longest period normal modes using recent large earthquakes, Geophys. Res. Lett.  38 L04303, doi:10.1029/2010GL046115. https://doi.org/10.1029/2010GL046115 Google Scholar CrossRef Search ADS   Deuss A., Ritsema J., van Heijst H., 2013. A new catalogue of normal-mode splitting function measurements up to 10 mHz, Geophys. J. Int.  193 920– 937. https://doi.org/10.1093/gji/ggt010 Google Scholar CrossRef Search ADS   Durek J.J., Romanowicz B., 1999. Inner core anisotropy inferred by direct inversion of normal mode spectra, Geophys. J. Int.  139 599– 622. https://doi.org/10.1046/j.1365-246x.1999.00961.x Google Scholar CrossRef Search ADS   Dziewonski A., Anderson D., 1981. Preliminary reference Earth model, Phys. Earth planet. Inter.  25 297– 356. https://doi.org/10.1016/0031-9201(81)90046-7 Google Scholar CrossRef Search ADS   Gharti H., Tromp J., 2017. A spectral-infinite-element solution of Poisson's equation: an application to self gravity , arXiv:1706.00855v1. Giardini D., Li X.-D., Woodhouse J., 1987. Three-dimensional structure of the Earth from splitting in free-oscillation spectra, Nature  325 405– 411. https://doi.org/10.1038/325405a0 Google Scholar CrossRef Search ADS   Giardini D., Li X.-D., Woodhouse J., 1988. Splitting functions of long-period normal modes of the Earth, J. geophys. Res.  93 13 716– 13 742. https://doi.org/10.1029/JB093iB11p13716 Google Scholar CrossRef Search ADS   Gilbert F., 1970. Excitation of the normal modes of the Earth by earthquake sources, Geophys. J. R. astr. Soc.  22 223– 226. https://doi.org/10.1111/j.1365-246X.1971.tb03593.x Google Scholar CrossRef Search ADS   Hara T., Tsuboi S., Geller R., 1991. Inversion for laterally heterogeneous earth structure using a laterally heterogeneous starting model: preliminary results, Geophys. J. Int.  104 523– 540. https://doi.org/10.1111/j.1365-246X.1991.tb05699.x Google Scholar CrossRef Search ADS   Hara T., Tsuboi S., Geller R., 1993. Inversion for laterally heterogeneous upper mantle S-wave velocity structure using iterative waveform inversion, Geophys. J. Int.  115 667– 698. https://doi.org/10.1111/j.1365-246X.1993.tb01487.x Google Scholar CrossRef Search ADS   He X., Tromp J., 1996. Normal-mode constraints on the structure of the Earth, J. geophys. Res.  101 20 053– 20 082. https://doi.org/10.1029/96JB01783 Google Scholar CrossRef Search ADS   Irving J., Deuss A., Andrews J., 2008. Wide-band coupling of Earth's normal modes due to anisotropic inner core structure, Geophys. J. Int.  174 919– 929. https://doi.org/10.1111/j.1365-246X.2008.03824.x Google Scholar CrossRef Search ADS   Ishii M., Tromp J., 1999. Normal-mode and free-air gravity constraints on lateral variations in velocity and density of Earth's mantle, Science  285 1231– 1236. https://doi.org/10.1126/science.285.5431.1231 Google Scholar CrossRef Search ADS PubMed  Karato S.-i., 1993. Importance of anelasticity in the interpretation of seismic tomography, Geophys. Res. Lett.  15 1623– 1626. https://doi.org/10.1029/93GL01767 Google Scholar CrossRef Search ADS   Koelemeijer P., Deuss A., Ritsema J., 2013. Observations of core-mantle boundary Stoneley modes, Geophys. Res. Lett.  40 2557– 2561. https://doi.org/10.1002/grl.50514 Google Scholar CrossRef Search ADS   Koelemeijer P., Deuss A., Ritsema J., 2017. Density structure of Earth's lowermost mantle from Stoneley mode splitting observations, Nat. Commun.  8 15241. https://doi.org/10.1038/ncomms15241 Google Scholar CrossRef Search ADS PubMed  Komatitsch D., Tromp J., 2002. Spectral-element simulations of global seismic wave propagation—I. Validation, Geophys. J. Int.  149 390– 412. https://doi.org/10.1046/j.1365-246X.2002.01653.x Google Scholar CrossRef Search ADS   Kuo C., Romanowicz B., 2002. On the resolution of density anomalies in the Earth's mantle using spectral fitting of normal-mode data, Geophys. J. Int.  150 162– 179. https://doi.org/10.1046/j.1365-246X.2002.01698.x Google Scholar CrossRef Search ADS   Lau H., Yang H.-Y., Tromp J., Mitrovica J., Latychev K., Al-Attar D., 2015. A normal mode treatment of semi-diurnal body tides on an aspherical, rotating and anelastic Earth, Geophys. J. Int.  202 1392– 1406. https://doi.org/10.1093/gji/ggv227 Google Scholar CrossRef Search ADS   Lau H., Mitrovica J., Davies J., Tromp J., Yang H.-Y., Al-Attar D., 2017. Tidal tomography constrains Earth's deep-mantle buoyancy, Nature  551 321– 326. https://doi.org/10.1038/nature24452 Google Scholar CrossRef Search ADS PubMed  Lay T., 2007. Deep earth structure—lower mantle and D″, in Treatise on Geophysics  Vol. 1, pp. 620– 654, eds Romanowicz B., Dziewoński A., Elsevier. Lognonné P., 1991. Normal modes and seismograms in an anelastic rotating Earth, J. geophys. Res.  96 20 309– 20 319. https://doi.org/10.1029/91JB00420 Google Scholar CrossRef Search ADS   Masters G., Laske G., Gilbert F., 2000. Matrix autoregressive analysis of free-oscillation coupling and splitting, Geophys. J. Int.  143 478– 489. https://doi.org/10.1046/j.1365-246X.2000.01261.x Google Scholar CrossRef Search ADS   Masters G., Barmine M., Kientz S., 2011. Mineos User Manual  Computational Infrastructure for Geodynamics, v. 1.0.2. Mosca I., Cobden L., Deuss A., Ritsema J., Trampert J., 2012. Seismic and mineralogical structures of the lower mantle from probabilistic tomography, J. geophys. Res.  117 B06304, doi:10.1029/2011JB008851. https://doi.org/10.1029/2011JB008851 Google Scholar CrossRef Search ADS   Moulik P., Ekström G., 2016. The relationships between large-scale variations in shear velocity, density, and compressional velocity in the Earth's mantle, J. geophys. Res.  121 2737– 2771. https://doi.org/10.1002/2015JB012679 Google Scholar CrossRef Search ADS   Pachhai S., Tkalčić H., Masters G., 2016. Estimation of splitting functions from Earth's normal mode spectra using the neighbourhood algorithm, Geophys. J. Int.  204 111– 126. https://doi.org/10.1093/gji/ggv414 Google Scholar CrossRef Search ADS   Park J., 1986. Synthetic seismograms from coupled free oscillations: Effects of lateral structure and rotation, J. geophys. Res.  91 6441– 6464. https://doi.org/10.1029/JB091iB06p06441 Google Scholar CrossRef Search ADS   Park J., 1987. Asymptotic coupled-mode expressions for multiplet amplitude anomalies and frequency shifts on an aspherical earth, Geophys. J. R. astr. Soc.  90 129– 169. https://doi.org/10.1111/j.1365-246X.1987.tb00679.x Google Scholar CrossRef Search ADS   Park J., 1990. The subspace projection method for constructing coupled-mode synthetic seismograms, Geophys. J. Int.  101 111– 123. https://doi.org/10.1111/j.1365-246X.1990.tb00761.x Google Scholar CrossRef Search ADS   Pekeris C., Jarosch H., 1958. The free oscillations of the earth, in Contributions in Geophysics in Honor of Beno Butenberg  eds Benioff H., Ewing M., Howell B., Press F., Pergamon Press. Resovsky J., Ritzwoller M., 1995. Constraining odd-degree Earth structure with coupled free oscillations, Geophys. Res. Lett.  22 2301– 2304. https://doi.org/10.1029/95GL01996 Google Scholar CrossRef Search ADS   Resovsky J., Ritzwoller M., 1998. New and refined constraints on three-dimensional Earth structure from normal modes below 3 mHz, J. geophys. Res.  103 783– 810. https://doi.org/10.1029/97JB02482 Google Scholar CrossRef Search ADS   Resovsky J., Ritzwoller M., 1999. Regularization uncertainty in density models estimated from normal mode data, Geophys. Res. Lett.  26 2319– 2322. https://doi.org/10.1029/1999GL900540 Google Scholar CrossRef Search ADS   Resovsky J., Trampert J., 2003. Using probabilistic seismic tomography to test mantle velocity-density relationships, Earth planet. Sci. Lett.  215, 121– 234. https://doi.org/10.1016/S0012-821X(03)00436-9 Google Scholar CrossRef Search ADS   Ritsema J., van Heijst H., Woodhouse J., 1999. Complex shear wave velocity structure imaged beneath Africa and Iceland, Science  286( 5446), 1925– 1931. https://doi.org/10.1126/science.286.5446.1925 Google Scholar CrossRef Search ADS PubMed  Ritsema J., van Heijst H., Woodhouse J., 2004. Global transition zone tomography, J. geophys. Res.  109 doi:10.1029/2003JB002610. https://doi.org/10.1029/2003JB002610 Ritsema J., Deuss A., van Heijst H., Woodhouse J., 2011. S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements, Geophys. J. Int.  184 1223– 1236. https://doi.org/10.1111/j.1365-246X.2010.04884.x Google Scholar CrossRef Search ADS   Rogister Y., Valette B., 2009. Influence of liquid core dynamics on rotational modes, Geophys. J. Int.  176 368– 388. https://doi.org/10.1111/j.1365-246X.2008.03996.x Google Scholar CrossRef Search ADS   Romanowicz B., 1987. Multiplet-multiplet coupling due to lateral heterogeneity: asymptotic effects on the amplitude and frequency of the Earth's normal modes, Geophys. J. R. astr. Soc.  90 75– 100. https://doi.org/10.1111/j.1365-246X.1987.tb00676.x Google Scholar CrossRef Search ADS   Tarantola A., Valette B., 1982. Generalized nonlinear inverse problems solved using the least squares criterion, Rev. Geophys. Space Phys.  20 219– 232. https://doi.org/10.1029/RG020i002p00219 Google Scholar CrossRef Search ADS   Trampert J., Deschamps F., Resovsky J., Yuen D., 2004. Probabilistic tomography maps chemical heterogeneities throughout the lower mantle, Science  306 853– 856. https://doi.org/10.1126/science.1101996 Google Scholar CrossRef Search ADS PubMed  Tromp J., Dahlen F., 1990. Free oscillations of a spherical anelastic earth, Geophys. J. Int.  103 707– 723. https://doi.org/10.1111/j.1365-246X.1990.tb05682.x Google Scholar CrossRef Search ADS   Um J., Dahlen F., 1992. Normal mode multiplet coupling on an aspherical, anelastic earth, Geophys. J. Int.  111 11– 34. https://doi.org/10.1111/j.1365-246X.1992.tb00551.x Google Scholar CrossRef Search ADS   Valentine A., Trampert J., 2012. Assessing the uncertainties on seismic source parameters: towards realistic error estimates for centroid-moment-tensor determinations, Phys. Earth. planet. Inter.  210–211 36– 49. https://doi.org/10.1016/j.pepi.2012.08.003 Google Scholar CrossRef Search ADS   Valentine A., Trampert J., 2016. The impact of approximations and arbitrary choices on geophysical images, Geophys. J. Int.  204 59– 73. Google Scholar CrossRef Search ADS   Woodhouse J., 1980. The coupling and attenuation of nearly resonant multiplets in the Earth's free oscillation spectrum, Geophys. J. R. astr. Soc.  61 261– 283. https://doi.org/10.1111/j.1365-246X.1980.tb04317.x Google Scholar CrossRef Search ADS   Woodhouse J., 1988. The calculation of the eigenfrequencies and eigenfunctions of the free oscillations of the earth and the sun, in Seismological Algorithms  pp. 321– 370, ed. Doornbos D., Academic Press. Woodhouse J., Dahlen F., 1978. The effect of a general aspherical perturbation on the free oscillations of the Earth, Geophys. J. R. astr. Soc.  53 335– 354. https://doi.org/10.1111/j.1365-246X.1978.tb03746.x Google Scholar CrossRef Search ADS   Woodhouse J., Deuss A., 2007. Earth's free oscillations, in Seismology and structure of the Earth  vol. 1 of Treatise on Geophysics, chap. 2, pp. 31– 65, eds Romanowicz B., Dziewoński A., Elsevier. Woodhouse J., Dziewonski A., 1984. Mapping the upper mantle: Three-dimensional modeling of Earth structure by inversion of seismic waveforms, J. geophys. Res.  89 5953– 5986. https://doi.org/10.1029/JB089iB07p05953 Google Scholar CrossRef Search ADS   Woodhouse J., Giardini D., 1985. Inversion for the splitting function of isolated low order normal mode multiplets, EOS, Trans. Am. geophys. Un.  66 300. Woodhouse J., Girnius T., 1982. Surface waves and free oscillations in a regionalized earth model, Geophys. J. R. astr. Soc.  68 653– 673. https://doi.org/10.1111/j.1365-246X.1982.tb04921.x Google Scholar CrossRef Search ADS   Woodhouse J., Giardini D., Li X.-D., 1986. Evidence for inner core anisotropy from free oscillations, Geophys. Res. Lett.  13 1549– 1552. https://doi.org/10.1029/GL013i013p01549 Google Scholar CrossRef Search ADS   Yang H.-Y., Tromp J., 2015. Synthetic free oscillation spectra: an appraisal of various mode-coupling methods, Geophys. J. Int.  203 1179– 1192. https://doi.org/10.1093/gji/ggv349 Google Scholar CrossRef Search ADS   APPENDIX: SYNTHETIC SEISMIC WAVEFORMS AND SPECTRA—THE NORMAL MODE FRAMEWORK A complete description of the theoretical basis for normal mode seismology lies well beyond the scope of this paper. Here, we aim to provide sufficient background to allow the reader to appreciate the motivations for—and distinctions between—self-, group- and full-coupling. In particular, we aim to highlight precisely where and why approximations need to be introduced within the various theories, providing insight into potential for systematic errors. Readers wishing for a more in-depth treatment are directed to a summary paper by Woodhouse & Deuss (2007), or the comprehensive account contained within Dahlen & Tromp (1998). Of course, our discussion here is very general: specific implementations of these techniques invariably adopt additional simplifications and restrictions, which may impact further upon accuracy. A1 The seismic wave equation The seismic wavefield within the Earth, s(x, t), can be described by the seismic wave equation. This has the general form   \begin{equation} \left(\mathcal {V}+\mathcal {W}\frac{\partial }{\partial t} +\rho \frac{\partial ^2}{\partial t^2}\right)\mathbf {s}=\mathbf {f} \end{equation} (A1)where $$\mathcal {V}$$ represents a linear integro-differential operator that depends upon the material properties of the Earth; $$\mathcal {W}$$ is an operator representing effects due to the Coriolis force; ρ = ρ(x) is the equilibrium density distribution within the Earth; and f = f(x, t) represents the force due to the seismic source. If viscoelastic effects are to be properly accounted for, the operator $$\mathcal {V}$$ is itself a function of time, and in this case $$\mathcal {V}\mathbf {s}$$ implies a convolution operation. Given specific realizations of the various quantities, we aim to solve eq. (A1) to obtain synthetic seismic observables. A general strategy for solving differential equations entails expanding the solution in terms of a known basis set, thereby reducing the equations to algebraic form. Although the solution remains the same irrespective of the basis chosen, the computational procedure required to obtain that solution can vary considerably. Thus, an appropriate choice of basis can be critical to developing efficient algorithms. In a seismological context, one strategy is to adopt a basis constructed on a finite-element–style mesh. This results in computational tools such as SPECFEM3D_GLOBE (e.g Komatitsch & Tromp 2002), where simulation costs scale with the length of time-series that must be calculated. Normal mode seismology takes a different approach: it relies upon basis functions that have global support, and the resulting algorithms have computational costs that depend largely on the desired frequency range of solution. Both approaches are formally complete and can provide faithful simulation of physical reality. However, we note that SPECFEM3D_GLOBE currently omits self-gravitation effects, as these are challenging to incorporate in a spectral element framework, and is therefore not suitable for modelling seismograms at the very low frequencies we consider in this paper (although efforts are being made to overcome this difficulty: see Gharti & Tromp 2017). A2 Spherical-earth eigenfunctions The particular family of basis functions we use are the eigenfunctions, or ‘normal modes’, of a spherically symmetric earth model. In this case, the material properties vary only with depth, and eq. (A1) simplifies considerably. The Coriolis term can be discarded, and the operator $$\mathcal {V}$$ takes a simpler form, which we denote by $$\bar{\mathcal {V}}$$. Similarly, we use $$\bar{\rho }$$ to highlight that the density distribution is also spherically symmetric. Neglecting the force term, and working in the frequency domain, it can then be shown that the differential equations admit a non-trivial solution, $$\boldsymbol {\mathcal {S}}_k(\mathbf {x})$$, at a discrete set of frequencies ωk. Thus, we have   \begin{equation} \left(\bar{\mathcal {V}}-\bar{\rho } \omega _k^2\right)\boldsymbol {\mathcal {S}}_k = 0 \end{equation} (A2)where $$\boldsymbol {\mathcal {S}}_k$$ represents the kth normal mode, and ωk is its corresponding eigenfrequency. This eigenvalue problem can be solved relatively straightforwardly (Pekeris & Jarosch 1958). There are an infinite number of $$(\omega _k, \boldsymbol {\mathcal {S}}_k)$$ pairs that satisfy eq. (A2), but only a finite set exist within any given frequency interval ω0 ≤ ωk ≤ ω1. Thus, it is feasible to compute a complete set of eigenfunctions within the given frequency range, and an algorithm described by Woodhouse (1988) provides a means for doing so. This method is implemented, for example, within the MINEOS software package described by Masters et al. (2011). It turns out (see e.g. Pekeris & Jarosch 1958) that distinct families of normal modes may be identified. Principally, two classes exist, known as ‘spheroidal’ and ‘toroidal’ modes; in addition, ‘radial’ modes are sometimes discussed as a separate category, although they are essentially a particular subset of the spheroidals. Furthermore, distinct sets of toroidal modes may be identified for each isolated solid region within the structural model, thus two groups of toroidal modes may be found for earth models such as PREM, where a solid inner core and mantle/crust region are separated by a fluid outer core. The modes take the spatial form of spherical harmonics laterally; their radial functions must be found by numerical integration. Thus, the index k on $${\boldsymbol {\mathcal {S}}}_k$$ may be recast as a set of four indices, k → (q, l, m, n). Here, q defines the ‘family’ of mode being considered; l and m are the angular degree and azimuthal order specifying the appropriate spherical harmonic Ylm(θ, ϕ); and n is the ‘overtone number’, which indexes the set of permissible radial functions. For present purposes, the main point of interest is that ωk depends only on q, l and n: in a spherically symmetric model, modes that differ only by their m-value share the same eigenfrequency. In the literature of normal mode seismology, it is common to see an individual (q, l, m, n) state described as a ‘singlet’; the (2l + 1) degenerate singlets corresponding to a given (q, l, n) triplet are then referred to as a ‘multiplet’. This triplet is often written in the form nql, with q being either S (in the case of spheroidal modes), or T (toroidals). For brevity, we will continue to make use of the index k within this discussion, unless there is reason to do otherwise. A2.1 Completeness and orthogonality The eigenfunctions of a spherically symmetric model can be shown to have some important properties, which stem from the self-adjointness of the operator $$\bar{ \mathcal {V}}$$. In particular, the infinite set of $${\boldsymbol {\mathcal {S}}}_k$$ form a complete basis for vector fields within the earth model, and they may be shown to be orthogonal, in the sense that   \begin{equation} \int _{V}\bar{\rho }(\mathbf {x}) \, {\boldsymbol {\mathcal {S}}}_j^*(\mathbf {x})\cdot {\boldsymbol {\mathcal {S}}}_k(\mathbf {x})\,\mathrm{d}^3\mathbf {x}=\delta _{jk} \end{equation} (A3)where δjk represents the Kronecker delta. Here, the integration volume V is the complete interior of the earth model, and an asterisk is used to denote complex conjugation. Strictly, these properties are only proven for entirely solid models, and we ignore any complications that might arise from the existence of ‘undertone modes’ in the Earth's fluid outer core. These are not thought to be problematic in a seismological context, as they exist at lower frequencies than we typically consider (e.g. Rogister & Valette 2009). A2.2 Seismograms in a spherically symmetric model Knowledge of the spherical-earth eigenfunctions allows straightforward synthesis of the wavefield expected within a spherically symmetric body. As shown by Gilbert (1970), and recounted in Woodhouse & Deuss (2007), the three-component wavefield can be reduced to the form   \begin{equation} {\bar{\mathbf {s}}}(\mathbf {x},t) = \sum _k \frac{1}{\omega _k^2}E_k (1-\exp \imath \omega _k t) \boldsymbol {\mathcal {S}}_k(\mathbf {x}) \end{equation} (A4)where the real part is understood and we only consider positive times. Ek represent the modal excitations, and may be obtained by calculating certain spatial integrals of the force term f(x, t). The process of computing seismograms by this formula is often referred to—for obvious reasons—as ‘mode summation’. In principle, the summation over k includes an infinite number of modes. However, we see that the dominant contribution of each mode to the seismogram lies at its eigenfrequency. Thus, if the seismogram is only required within a finite frequency band, only modes lying within that band need to be summed, and this makes the summation finite. In practice, we always work with band-limited seismic data, which is a necessary consequence of digital signal recording. If the maximum frequency present in our data set is ωmax, it is straightforward—and sufficient—to compute a mode catalogue containing all modes with ωk ≤ ωmax. A3 Seismograms in an aspherical earth model We now return to consider general, aspherical earth models, and the solution of eq. (A1). Invariably, the aspherical model will be specified as a 3-D structure superimposed upon a spherically symmetric reference model (often based on the model PREM, as described in Dziewonski & Anderson 1981). In what follows, we assume that the $$\boldsymbol {\mathcal {S}}_k$$ have been obtained using the relevant spherically symmetric reference structure. In addition, we implicitly assume that all interfaces within the aspherical model remain concentric spheres: in other words, no topography is present, either on the free surface or on internal discontinuities (such as the core–mantle boundary). To incorporate such topography, studies have traditionally relied on first-order boundary perturbation theory (Woodhouse & Dahlen 1978; Woodhouse 1980), which is inherently an approximation. Recently, an exact treatment of boundary perturbations has been developed by Al-Attar & Crawford (2016), although this theory has yet to be fully implemented. In any case, the effects of realistic-scale topography are not likely to impact the conclusions of this paper, and need not be considered further here. Given that the $$\boldsymbol {\mathcal {S}}_k$$ form a complete basis, we can express the wavefield s(x, t) in the form   \begin{equation} \mathbf {s}(\mathbf {x},t) = \sum _k u_k(t) {\boldsymbol {\mathcal {S}}}_k(\mathbf {x})\mathrm{,} \end{equation} (A5)where the real part is again understood. uk represent time-varying expansion coefficients. Strictly, this is only valid if the summation includes the entire and infinite set of eigenfunctions for the spherically symmetric model. However, we note the similarity with eq. (A4), and the fact that realistic aspherical models contain only relatively weak perturbations to the spherically symmetric reference model. It therefore seems reasonable to assert that modes with eigenfrequency far in excess of ωmax are unlikely to contribute within the frequency band of interest. We introduce a cut-off frequency, ωc ≥ ωmax, and work with the finite set of $$\boldsymbol {\mathcal {S}}_k$$ having ωk ≤ ωc. A3.1 Full-coupling In order to make use of this expansion, eq. (A5), we need to obtain expressions for the coefficients uk. Substituting into the wave equation, eq. (A1), we find   \begin{equation} \sum _k \left(u_k\mathcal {V}{\boldsymbol {\mathcal {S}}}_k+\dot{u}_k\mathcal {W}\boldsymbol {\mathcal {S}}_k+\rho \ddot{u}_k{\boldsymbol {\mathcal {S}}}_k \right)= \mathbf {f} \end{equation} (A6)and therefore, for any eigenfunction $$\boldsymbol {\mathcal {S}}_j$$, we obtain   \begin{eqnarray} &&{\sum _k\left(u_k \int _V \boldsymbol {\mathcal {S}}^*_j\cdot \mathcal {V}\boldsymbol {\mathcal {S}}_k\,\mathrm{d}^3\mathbf {x} + \dot{u}_k\int _{V}\boldsymbol {\mathcal {S}}^*_j\cdot \mathcal {W}\boldsymbol {\mathcal {S}}_k\,\mathrm{d}^3 \mathbf {x}\right.}\nonumber\\ &&{\quad+\,\,\left.\ddot{u}_k\int _{V}\rho \,\boldsymbol {\mathcal {S}}^*_j\cdot \boldsymbol {\mathcal {S}}_k\,\mathrm{d}^3\mathbf {x}\right) = \int _{V}\boldsymbol {\mathcal {S}}_j^* \cdot \mathbf {f}\,\mathrm{d}^3 \mathbf {x}.} \end{eqnarray} (A7)Since the sum only contains a finite number of terms corresponding to frequencies lower than ωc, this can equivalently be expressed as a matrix equation,   \begin{equation} \mathbf {Vu} + \mathbf {W\dot{u}}+\mathbf {P\ddot{u}}= \mathbf {q} \end{equation} (A8)where the elements of P, W, V and q are obtained by evaluating the appropriate integrals from eq. (A7). Taking the Fourier transform allows us to eliminate the derivatives, and introduce the matrix M(ω), such that   \begin{equation} \left[\mathbf {\tilde{V}}(\omega ) +\imath \omega \mathbf {W}-\omega ^2 \mathbf {P}\right]\mathbf {\tilde{u}} = \mathbf {M}(\omega ) \mathbf {\tilde{u}} = \mathbf {\tilde{q}} \end{equation} (A9)where the frequency-dependence of $$\mathbf {\tilde{V}}$$ arises from the inclusion of viscoelasticity within $$\mathcal {V}$$; strictly, eq. (A8) involves a convolution for the first term. The elements of M and $$\mathbf {\tilde{q}}$$ are readily evaluated, and thus obtaining $$\mathbf {\tilde{u}}$$ is a straightforward exercise in linear algebra. The ‘full-coupling’ approach involves solving eq. (A9) exactly. In reaching this point, we have introduced only one approximation: we have truncated the modal expansion at frequency ωc. To date, there has been little effort to investigate the errors introduced by this truncation, or assess how ωc should be chosen. Most of the existing studies that employ full-coupling assume it is sufficient to use ωc = ωmax, but there is little experimental or theoretical justification to support this. Thus, one goal of this paper is to explore how changes to ωc impact the seismograms produced, and to identify where the cut-off must be placed to ensure computations are essentially exact at frequencies below ωmax. We note that a number of ‘full-coupling’ studies introduce additional approximations within eq. (A9)—similar to those discussed in Appendix A4, below—designed to reduce computational complexity. It has been shown that these introduce only minor errors in results (Deuss & Woodhouse 2001). However, it is possible to avoid these simplifications, and the full-coupling calculations performed in this paper are based on solving eq. (A9) exactly. To achieve this, we adopt an iterative, preconditioned method described in Al-Attar et al. (2012), building on earlier work by Hara et al. (1993) and Deuss & Woodhouse (2004). A3.2 Self- and group-coupling Although solving eq. (A9) is conceptually straightforward, computational costs grow rapidly as ωc increases. Broadly speaking, the number of modes having ωk ≤ ωc (and thus, the dimension of the system in eq. (A9)) scales proportionally with $$\omega _c^3$$; the cost of obtaining $$\mathbf {\tilde{u}}$$ then grows as $$\omega _c^6$$. Computational costs rapidly become prohibitive, even with modern computational resources; for early studies, full-coupling was intractable. In consequence, it was necessary to identify strategies for simplifying the computational task. Looking again at eq. (A7), we can consider how the equations behave in the absence of any aspherical structure. In this case, all the integrals vanish unless j = k and thus, in the limit of a purely spherical model, the matrix M has non-zero elements only along its diagonal. This leads directly to our expression for seismograms in a spherically symmetric earth, eq. (A4), as is expected. However, it also motivates the assumption—borne out by empirical evidence—that for typical weakly aspherical models, M remains diagonally dominant. Elements lying far from the diagonal are, in general, relatively small—and in group- and self-coupling, they are therefore ignored. In eq. (A9), a general element of M is defined   \begin{equation} M_{jk}(\omega ) = \int _{V} \boldsymbol {\mathcal {S}}_j^*(\mathbf {x})\cdot \left[\mathcal {V} + \imath \omega \mathcal {W} - \omega ^2\rho (\mathbf {x})\right]\boldsymbol {\mathcal {S}}_k(\mathbf {x})\,\mathrm{d}^3 \mathbf {x}. \end{equation} (10a)In self-coupling, this is modified by introducing the assertion that   \begin{equation} M_{jk} = M_{(q_j,l_j,n_j)(q_k,l_k,n_k)}= 0 \,\,\textrm{unless}\,\, \left\lbrace \begin{array}{@{}l@{\quad }l@{}}q_j=q_k\\ l_j=l_k\\ n_j=n_k.\end{array}\right. \end{equation} (10b)In other words, self-coupling computes integrals using only the constituent singlets of each multiplet, and assumes that interactions between different multiplets can be entirely ignored. In group-coupling, the multiplets are assigned to groups. We introduce κ to denote a given multiplet, κ = {(q, l, m, n): m ∈ [−l, l]}. We then introduce groups, $$\mathcal {G}$$, containing one or more multiplets: $$\mathcal {G}_i = \lbrace \kappa _1^{(i)}, \kappa _2^{(i)},\ldots \rbrace$$. In general, all multiplets within a group will have similar eigenfrequency, although individual authors may differ over the particular groupings adopted. Group-coupling then replaces eq. (A10b) by the assertion that   \begin{equation} M_{jk} = 0 \,\,\textrm{unless}\,\, \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(q_j,l_j,m_j,n_j)\in \mathcal {G}\\ (q_k,l_k,m_k,n_k)\in \mathcal {G} \end{array}\right. \end{equation} (10c) for a single group, $$\mathcal {G}$$. Thus, some limited interactions are allowed between multiplets. However, (assuming a sensible ordering of the modes), both self- and group-coupling cause the matrix M to take on a block-diagonal form. Exploiting this structure allows a great reduction in the computational effort required to solve eq. (A9), at the expense, of course, of a deterioration in accuracy. A4 Normal modes in aspherical models Thus far, we have discussed the use of full-, group- and self-coupling to generate seismic spectra in aspherical earth models. These techniques make use of the eigenfunctions and eigenfrequencies of spherically symmetric models. However, it is also possible to compute the normal modes of general, aspherical models. One route to doing so is to adopt the same general strategy as we have already outlined: assume that the aspherical eigenfunction can be expressed in terms of the modes of a spherically symmetric reference model, and then solve for the necessary coefficients. This results in essentially the same analysis as in Appendix A3.1, except with time-invariant uk, and no force term. One may then regard eq. (A9) as a nonlinear eigenvalue problem, which can (in principle) be solved to yield the necessary expansion coefficients. The procedure is computationally challenging, limiting its range of application, especially in early studies. To simplify matters, we can make the same assumption as in self-coupling: that integrals of the form $$\int _V \boldsymbol {\mathcal {S}}_i^*(\mathbf {x})\cdot \mathcal {A}(\omega )\boldsymbol {\mathcal {S}}_j(\mathbf {x})\,\mathrm{d}^3 \mathbf {x}$$ are only non-zero when the singlets $$\boldsymbol {\mathcal {S}}_i$$ and $$\boldsymbol {\mathcal {S}}_j$$ correspond to the same multiplet. This allows the eigenvalue problem of eq. (A9) to be decomposed into a sequence of independent calculations, one for each multiplet. Furthermore, a perturbation-theory approach may be used to approximate the nonlinear eigenvalue problem as one in standard form. Given the orthogonality relationship, eq. (A3), we can write P = I + P1, where I is an identity matrix and P1 is the part arising only from the aspherical density structure (i.e. $$\rho -\bar{\rho }$$). We then assume that M(ω) can be expanded around a reference eigenfrequency of the multiplet, ω0, such that (cf. eq. A9)   \begin{equation} \mathbf {M}(\omega )\approx \mathbf {\tilde{V}}^{(\kappa \kappa )}(\omega _0)+\imath \omega _0\mathbf {W}^{(\kappa \kappa )} -\omega _0^2 \mathbf {P_1}^{(\kappa \kappa )} - \omega ^2\mathbf {I} \end{equation} (A11)where the superscript (κκ) is used to emphasize that each multiplet is being treated independently. This approximation relies on the assumption that |ω − ω0| is small, and that the earth model is dominated by spherical structure, allowing various additional terms involving P1, W and $$\partial \mathbf {\tilde{V}}/\partial \omega$$ can be neglected. It is common to introduce a ‘splitting matrix’, $$\mathbf {H} = \mathbf {\tilde{V}}(\omega _0)+\imath \omega _0\mathbf {W} -\omega _0^2 \mathbf {P_1}$$, allowing eq. (A9) to be written in the form   \begin{equation} \left(\mathbf {H}^{(\kappa \kappa )} - \omega ^2 \mathbf {I}\right)\mathbf {\tilde{u}}^{(\kappa \kappa )} = \mathbf {\tilde{q}}. \end{equation} (A12)Thus, the contribution of each multiplet to spectra—within the self-coupling approximation—is then found by diagonalizing the (2l + 1) × (2l + 1) matrix H. The (2l + 1) constituent singlets within the multiplet each contribute at a slightly different eigenfrequency, removing the degeneracy seen in spherically symmetric models. Thus, the spectral contribution of each multiplet becomes a set of (2l + 1) closely spaced—but distinct—peaks. This phenomenon is known as ‘splitting’. In practice, the resolution of spectral measurements may not allow individual singlets to be distinguished, but the overall shape of the peak corresponding to each multiplet depends upon the splitting, and thus upon the detailed properties of the Earth model. In the time domain, the contribution of the multiplet to seismograms can be written in the form   \begin{equation} \mathbf {s}^{(\kappa \kappa )}(\mathbf {x},t) = \mathbf {r}(\mathbf {x})\cdot \exp \left(\imath \mathbf {H}^{(\kappa \kappa )}t\right)\cdot \mathbf {v} \end{equation} (A13)where v and r are vectors that may be computed for the source and receiver, respectively; this result is due to Woodhouse & Girnius (1982). The discussion in this section has focused on self-coupling involving only individual isolated multiplets. Similar results can be obtained for group-coupling–assumptions: in effect, group-coupling can be regarded as self-coupling of a ‘super-multiplet’, allowing our analysis to translate straightforwardly. However, defining the reference frequency ω0 becomes more challenging, as the various multiplets within the group have different eigenfrequencies. Typically, some sort of average of these eigenfrequencies is used (for a recent discussion on specific averages, see Deuss & Woodhouse 2001). We note that the requirement for |ω − ω0| to remain small acts to restrict the size of the coupling group: only multiplets with similar eigenfrequency may be used. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Exact free oscillation spectra, splitting functions and the resolvability of Earth's density structure

Loading next page...
 
/lp/ou_press/exact-free-oscillation-spectra-splitting-functions-and-the-hN0u9zr3yj
Publisher
The Royal Astronomical Society
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx539
Publisher site
See Article on Publisher Site

Abstract

Summary Seismic free oscillations, or normal modes, provide a convenient tool to calculate low-frequency seismograms in heterogeneous Earth models. A procedure called ‘full mode coupling’ allows the seismic response of the Earth to be computed. However, in order to be theoretically exact, such calculations must involve an infinite set of modes. In practice, only a finite subset of modes can be used, introducing an error into the seismograms. By systematically increasing the number of modes beyond the highest frequency of interest in the seismograms, we investigate the convergence of full-coupling calculations. As a rule-of-thumb, it is necessary to couple modes 1–2 mHz above the highest frequency of interest, although results depend upon the details of the Earth model. This is significantly higher than has previously been assumed. Observations of free oscillations also provide important constraints on the heterogeneous structure of the Earth. Historically, this inference problem has been addressed by the measurement and interpretation of splitting functions. These can be seen as secondary data extracted from low frequency seismograms. The measurement step necessitates the calculation of synthetic seismograms, but current implementations rely on approximations referred to as self- or group-coupling and do not use fully accurate seismograms. We therefore also investigate whether a systematic error might be present in currently published splitting functions. We find no evidence for any systematic bias, but published uncertainties must be doubled to properly account for the errors due to theoretical omissions and regularization in the measurement process. Correspondingly, uncertainties in results derived from splitting functions must also be increased. As is well known, density has only a weak signal in low-frequency seismograms. Our results suggest this signal is of similar scale to the true uncertainties associated with currently published splitting functions. Thus, it seems that great care must be taken in any attempt to robustly infer details of Earth's density structure using current splitting functions. Computational seismology, Seismic tomography, Surface waves and free oscillations, Theoretical seismology 1 INTRODUCTION Our understanding of the Earth's large-scale interior structure and dynamics draws heavily on observations of seismic free oscillations (‘normal modes’). In particular, measurements of free-oscillation spectra at long periods represent one of the few classes of geophysical observable with meaningful sensitivity to lateral variations in density within the Earth. As such, free oscillation spectra have played a central role in attempts to determine the nature of the ‘Large Low Shear Velocity Provinces’ that have been identified within the lowermost mantle, and in attempts to identify the relative importance of thermal and chemical effects as driving forces for mantle convection (e.g. Ishii & Tromp 1999; Trampert et al.2004; Lay 2007; Mosca et al.2012; Moulik & Ekström 2016; Koelemeijer et al.2017). However, as this paper will show, computational limitations have led most such studies to rely upon seismological approximations that might be inadequate for imaging density variations. As such, their conclusions must be approached with considerable caution. In order to derive information about Earth's interior structure from seismic data, we must search for models which yield synthetic (predicted) observables which agree with those actually measured. Many strategies exist for performing this search, but all rely on the presumption that synthetic observables are computed accurately. The theoretical basis for computing synthetic normal mode spectra is well-developed, relying on a technique known as ‘normal mode coupling theory’ (e.g. Dahlen 1968, 1969; Woodhouse & Dahlen 1978; Woodhouse 1980; Woodhouse & Girnius 1982; Park 1986, 1987, 1990; Romanowicz 1987; Tromp & Dahlen 1990; Lognonné 1991; Hara et al.1991, 1993; Um & Dahlen 1992; Deuss & Woodhouse 2001; Al-Attar et al.2012; Yang & Tromp 2015; Lau et al.2015). Unlike some earlier studies, our calculations are based on the work of Al-Attar et al. (2012) and take the full effect of 3-D density variations into account, following the theory developed in Woodhouse (1980). Our calculations also allow us to take attenuation fully into account, although in the examples shown in this paper, only the spherically symmetric attenuation profile of PREM (Dziewonski & Anderson 1981) is used. Thus, our framework is physically complete, and allows the seismic response of any Earth-like body to be expressed exactly. However, the resulting expressions involve an infinite series expansion. This must be truncated in any computational implementation of the theory, introducing a truncation error within the synthetic spectra. Various common truncation schemes exist, motivated from physical principles. Computational costs scale unpleasantly with the number of terms retained in the series, so early studies had little choice but to adopt simplifying strategies. The most limiting approximation is known as ‘self-coupling’, while ‘group-coupling’ incorporates additional terms (we explain the distinction in more detail below). However, a number of studies, starting with Deuss & Woodhouse (2001), have demonstrated that neither can be relied upon within the frequency range of interest. In particular, Al-Attar et al. (2012) highlighted that use of either the self- or group-coupling approximations when computing synthetic spectra leads to errors of similar magnitude to the signal likely to be attributable to lateral variations in density within the Earth. Subsequently, a comprehensive study by Yang & Tromp (2015) concluded that self-coupling is ‘marginally acceptable’ when considering spectra below around 1.5 mHz, and ‘unacceptable’ above this point. Group-coupling was found to also yield significant errors, albeit at a lower level than in self-coupling. In principle, one could envisage improving group-coupling strategies by computing coupling strengths between modes—perhaps based on the work of Park (1987)—and ensuring that all significant interactions are accounted for. However, the notion of coupling strength is itself based on single-scattering approximations, and it is not clear how accurate the results would be. The optimal approach is therefore to adopt a strategy known as ‘full-coupling’. Despite the name, this continues to involve a truncation of the infinite series, but whereas self- and group-coupling are based around the presumption that the series is dominated by only a few terms, full-coupling aims to retain all terms that have any potential to contribute significantly; again, a more precise definition will be given in due course. Building on the existing body of work, the first aim of the present paper is to identify more precisely the conditions under which full-coupling can be regarded as ‘sufficiently accurate’. This information then permits us to perform high-quality synthetic experiments to address the second question: can robust information about earth structure be inferred from measurements of ‘splitting functions’ (Woodhouse & Giardini 1985; Giardini et al.1987, 1988)? These are a particular representation of information derived from seismic spectra, and they are an important ingredient in many current models of Earth structure (e.g. Ritsema et al.1999; Ishii & Tromp 1999; Trampert et al.2004; Ritsema et al.2011; Mosca et al.2012; Moulik & Ekström 2016; Koelemeijer et al.2017). However, the theory that connects splitting functions to structural parameters relies upon precisely the same assumptions as are inherent to self- and group-coupling. Given that these have been demonstrated to introduce significant bias into synthetic spectra, are splitting function inferences meaningful? The answer is not a foregone conclusion—the validity of any approximation is entirely context-specific—but the question must clearly be addressed. The first to address this question were Resovsky & Ritzwoller (1998). Being limited in computational resources, they estimated the theoretical error by a transfer function. Having at our disposal an arbitrarily precise forward solver, we can, for the first time, quantify this theoretical error properly. We therefore perform a simple, but definitive test: we compute accurate synthetic seismograms for the earth model S20RTS (Ritsema et al.1999), making use of full-coupling. We then measure splitting functions from this data set as if it were observed data. We also compute the ‘true’ splitting functions, by calculating the appropriate radial integrals for the known input model. Comparing the two, we find noticeable discrepancies. Repeating the experiment for synthetic spectra using the self-coupling approximation allows us to separate the contribution from theoretical approximations and regularization in the uncertainty of the measured splitting functions. Finally repeating the experiment for models where density perturbations are set to zero, we conclude that splitting functions are contaminated by a significant approximation error, a similar regularization error and contain a density signal of the level of the total measurement uncertainty. 2 THEORETICAL BACKGROUND Imaging the Earth's interior involves searching for the model parameters which can match observed seismic data. We therefore require a forward model, or a way to compute the observations that we would expect to see if the earth had any given structure. Ideally, this forward model needs to encapsulate the full physics of wave propagation. If it does not, the seismograms predicted for each candidate structure will be inaccurate, and the results of the imaging process may be biased by this systematic error (for a full discussion, see Valentine & Trampert 2016). 2.1 The forward problem ‘Normal mode seismology’ represents one framework for constructing a seismological forward model, particularly suited to low-frequency, global-scale simulations. The theory incorporates all important physical effects, and can be used to produce seismograms that faithfully reproduce physical reality. However, the computational resources required are significant, rendering full calculations intractable for early studies. As a result, a variety of simplifications were developed. These can greatly reduce the computational costs, at the expense of a degradation in accuracy. A brief, self-contained account of the underlying theory is given in Appendix, emphasizing how and why these various approximations can be introduced. Here, we sketch some key points. 2.1.1 Reference eigenfunctions For any spherically symmetric reference model (e.g. PREM, Dziewonski & Anderson 1981), it is relatively straightforward to compute eigenfunctions $$\boldsymbol {\mathcal {S}}_k$$ and their associated eigenfrequencies ωk; each pair $$(\boldsymbol {\mathcal {S}}_k,\omega _k)$$ represents one ‘normal mode’ (e.g. Woodhouse 1988). Conventionally, four indices are used to describe the modes: k → (q, l, m, n), with q denoting the ‘class’ (principally ‘spheroidal’, q = S, or ‘toroidal’, q = T ); l, the angular degree, and m, the azimuthal order, identifying the spherical harmonic Ylm(θ, ϕ) which describes the lateral form of $$\boldsymbol {\mathcal {S}}_k$$; and n, the overtone number, identifying its radial behaviour. The eigenfrequency can be shown to depend only upon q, l and n. For each (q, l, n) combination, there are therefore 2l + 1 eigenfunctions that have a different spatial form, but share a common eigenfrequency. This degenerate set is often referred to as a ‘multiplet’, and identified by a label in the form nql(e.g. 0S2). The constituent eigenfunctions of each multiplet are then referred to as ‘singlets’. 2.1.2 Computing synthetic seismograms A central theme of normal mode seismology is that these eigenfunctions, obtained for a spherically symmetric reference model, may be used as a convenient set of basis functions for representing vector fields within the (laterally heterogeneous) Earth. In particular, the seismic waveforms s(x, t) can be expressed in the form (cf. eq. A5)   \begin{equation} \mathbf {s}(\mathbf {x},t) = \sum _k u_k(t)\boldsymbol {\mathcal {S}}_k(\mathbf {x}) \end{equation} (1)where $$\boldsymbol {\mathcal {S}}_k$$ is a vector field representing an eigenfunction of the reference model. There are an infinite number of eigenfunctions $$\boldsymbol {\mathcal {S}}_k$$, and so in principle this representation of the seismic wavefield requires an infinite summation. However, if seismograms are to be filtered so that they contain no component above a given frequency, ωmax, it is justifiable to neglect modes with eigenfrequency ωk ≫ ωmax. To make this concrete, we define a cut-off frequency, ωc ≥ ωmax, and work with the finite set of modes having ωk ≤ ωc. Then, the Fourier-transformed expansion coefficients $$\tilde{u}_k(\omega )$$ may be found by solving a system of linear equations (cf. eq. A9)   \begin{equation} \mathbf {M}(\omega ) \mathbf {\tilde{u}} = \mathbf {\tilde{q}} \end{equation} (2)Here, M(ω) is a frequency-dependent ‘coupling matrix’ which can be computed for any earth model presenting 3-D variations in density, elastic constants and attenuation, while the vector $$\mathbf {\tilde{q}}$$ is a representation of the seismic source (typically an earthquake). 2.1.3 Full-coupling The term ‘full-coupling’ implies that we solve this linear system (i.e. eq. 2) directly. If we do this, we obtain seismograms that are essentially physically complete. They incorporate only one fundamental approximation: truncation of the infinite series at frequency ωc. The first goal of this paper is to quantify the effect of this truncation upon synthetic seismic spectra, and to assess where ωc should lie in order for the truncation to have no appreciable impact upon calculations for frequencies below ωmax. 2.1.4 The self- and group-coupling approximations As ωc increases, the dimension of the matrix M increases, and thus the cost of solving eq. (2) grows rapidly. The elements of M are found by computing integrals of the form (cf. eq. A10)   \begin{equation} M_{ij}(\omega ) = \int _V \boldsymbol {\mathcal {S}}_i^*(\mathbf {x})\cdot \mathcal {A}(\omega )\boldsymbol {\mathcal {S}}_j(\mathbf {x})\,\mathrm{d}^3\mathbf {x} \end{equation} (3a)where $$\mathcal {A}(\omega )$$ is an integro-differential operator that depends upon frequency, and also upon the earth model. The elements tend to be relatively small unless $$\boldsymbol {\mathcal {S}}_i$$ and $$\boldsymbol {\mathcal {S}}_j$$ have similar eigenfrequencies. This leads to structure within M, linked to the degeneracy within each multiplet; it is dominated by the elements lying close to the diagonal. Motivated by this, the ‘self-coupling approximation’ modifies the definition of M, introducing an over-riding assertion that   \begin{equation} M_{ij}(\omega ) = 0 \,\,\textrm{if}\,\,\boldsymbol {\mathcal {S}}_i\ \textrm{and}\ \boldsymbol {\mathcal {S}}_j\,\,{\rm come\ from\ different\ multiplets.} \end{equation} (3b)This results in M taking a block-diagonal form, with each multiplet nql contributing a block of dimension (2l + 1) × (2l + 1). Solution of eq. (2) is then much easier, as it decouples into separate calculations for each multiplet. Of course, the resulting solution is only an approximation to the one that would be obtained if the full-coupling matrix were used. The ‘group-coupling’ approximation allows for limited interaction between multiplets that are adjacent in frequency. The spectrum is divided into ‘groups’, with each group containing one or more multiplets. Studies may differ in the groupings used. Then, eq. (3b) is relaxed slightly, so that instead   \begin{equation} M_{ij}(\omega )= 0 \,\,\textrm{unless}\,\,\boldsymbol {\mathcal {S}}_i\ {\rm and}\ \boldsymbol {\mathcal {S}}_j\,{\rm belong\, to\, the\, same\, group.} \end{equation} (3c) Again, this simplifies the solution of eq. (2), allowing separate computations for each group. Studies have shown that the resulting seismograms are more accurate than those of self-coupling, but errors remain significant (e.g. Deuss & Woodhouse 2001; Al-Attar et al.2012; Yang & Tromp 2015). 2.1.5 Normal modes of an aspherical earth model So far, we have outlined how to obtain seismograms in heterogeneous earth models. It is also possible to compute the normal modes of a general (non-spherically symmetric) earth model (see Appendix A4) and again, group- or self-coupling approximations may be adopted if required. We find that the modes are no longer degenerate: each singlet within a multiplet has a distinct eigenfrequency. This property is sometimes referred to as ‘splitting’ of the multiplet, relative to the spherically symmetric reference. Given that realistic earth models contain only weak lateral heterogeneity, the typical frequency shifts involved are small, and due to finite-length time series, we do not generally expect to be able to resolve isolated singlets in spectral observations. Nevertheless, aspherical structure manifests itself in the overall shape of each ‘peak’ in the spectrum of a seismogram. 2.2 From spectra to structure: the splitting functions Given the ability to compute synthetic spectra, one may contemplate formulating an inverse problem that aims to infer earth structure from observations. This could be implemented in a variety of ways, each with pros and cons. One popular approach involves the use of ‘splitting functions’. These were introduced within Woodhouse & Giardini (1985), and developed in detail by Giardini et al. (1987, 1988). These original papers are developed under the self-coupling assumption, where each multiplet is regarded as entirely isolated from every other. This was recognized to be a simplification of the underlying physics, and therefore an approximation, but it was necessary if calculations were to be feasible using the computational resources then available. Making this approximation, the contribution of a given multiplet κ = nql within seismograms is entirely characterized by the (2l + 1) × (2l + 1) diagonal sub-block in the coupling matrix arising from eqs. (3a & 3b), which we denote M(κκ). Similarly, the multiplet can be directly identified with a specific portion of an observed seismic spectrum. Thus, given a suite of observations of the relevant portion at different locations and from various seismic events, one can formulate an inverse problem to estimate the elements of M(κκ). Splitting functions provide a mechanism for parametrizing this inversion, designed to provide a straightforward link between M(κκ) and the underlying earth model. 2.2.1 Splitting functions for isolated groups of multiplets The splitting function framework was extended by Resovsky & Ritzwoller (1995) and others, allowing for cross-coupling between specific pairs (or small groups) of multiplets with similar frequencies as in group-coupling. For simplicity, we describe here the case where two multiplets, κ = nql and $$\kappa ^\prime = {}_{n^\prime }q^\prime _{l^\prime }$$, are allowed to interact with one another, but are assumed to be isolated from the remainder of the spectrum. Extension to larger groups, or restriction to the self-coupling case, follows straightforwardly. The isolated-group assumption implies that the only relevant part of M(ω) is the symmetric diagonal block   \begin{equation} \left(\begin{array}{@{}cc@{}}\mathbf {{M}}^{(\kappa \kappa )} & \mathbf {{M}}^{(\kappa \kappa ^\prime )}\\ {\mathbf {{M}}^{(\kappa ^\prime \kappa )}}&\mathbf {{M}}^{(\kappa ^\prime \kappa ^\prime )}\end{array}\right). \end{equation} (4)As discussed in Appendix A4, some additional assumptions allow us to introduce the splitting matrix, H, such that $$\mathbf {{M}}^{(\kappa \kappa ^\prime )} \approx \mathbf {H}^{(\kappa \kappa ^\prime )}-\omega ^2\mathbf {I}$$. We then parametrize this by introducing a set of complex-valued ‘splitting coefficients’ cst such that   \begin{equation} \mathbf {H}^{(\kappa \kappa ^\prime )}_{mm^\prime } = \omega _0^2\delta _{\kappa \kappa ^\prime } + \omega _0 \mathbf {W}_{mm^\prime }^{(\kappa \kappa ^\prime )}+\sum _{s=l - l^\prime }^{l+l^\prime }\sum _{t=-s}^s\gamma _{ll^\prime s}^{mm^\prime t} c_{st}^{(\kappa \kappa ^\prime )} \end{equation} (5)where the γ are defined in terms of spherical harmonic triple products,   \begin{equation} \gamma _{ll^\prime s}^{mm^\prime t} = \int _0^{2\pi }\int _0^\pi Y^*_{lm}(\theta ,\phi )Y_{l^\prime m^\prime }(\theta ,\phi )Y_{st}(\theta ,\phi )\sin \theta \,\mathrm{d}\theta \mathrm{d}\phi . \end{equation} (6)This integral can be expressed straightforwardly in terms of the Wigner-3j symbols, and vanishes for many combinations of spherical harmonics. Together, the coefficients describe a ‘splitting function’   \begin{equation} \eta ^{(\kappa \kappa ^\prime )}(\theta ,\phi ) = \sum _{s=l-l^\prime }^{l+l^\prime }\sum _{t=-s}^s c_{st}^{(\kappa \kappa ^\prime )} Y_{st}(\theta ,\phi ) \end{equation} (7)which may be regarded as a map depicting spatial variations in spectral properties. Given the symmetry of M, three splitting functions are required to fully characterize our isolated pair of multiplets: η(κκ), $$\eta ^{(\kappa ^\prime \kappa ^\prime )}$$ and $$\eta ^{(\kappa \kappa ^\prime )}$$. Following Woodhouse & Girnius (1982), the contribution of the two multiplets within each seismogram can be computed via diagonalization of eq. (4) (see eq. A13 and the associated discussion). Thus, it is possible to pose the inverse problem whereby the cst are determined from observed data; typically, the fiducial frequency ω0 and parameters describing attenuation (which we will not discuss further) are also refined within this process. The inversion is found to be weakly non-linear, and is typically tackled using a Bayesian formulation of the least-squares algorithm (as per Tarantola & Valette 1982). Invariably, the data is insufficient to uniquely constrain a solution, and it is therefore necessary to introduce prior information (i.e. regularization). Discussion of the extent to which this influences results may be found in numerous sources (e.g. Resovsky & Ritzwoller 1999; Kuo & Romanowicz 2002; Pachhai et al.2016). In the isolated-group approximation, the splitting coefficients can be straightforwardly related to structural parameters. For the purposes of illustration, we assume that the aspherical earth model specifies P-wave speed (α), S-wave speed (β), and density (ρ) perturbations at every point with respect to a spherically symmetric reference model. Furthermore, we assume that the model is laterally parametrized in terms of spherical harmonics to maximum degree L. Thus, the P-wave speed field may be expressed   \begin{equation} \alpha (r,\theta ,\phi ) = \sum _{l=0}^L \alpha _{lm}(r) Y_{lm}(\theta ,\phi ) \end{equation} (8)with identical expressions for β and ρ. By deriving the exact expression for $$\mathbf {H}^{(\kappa \kappa ^\prime )}_{mm^\prime }$$ (e.g. Woodhouse 1980), and then comparing this with eq. (5), it is possible to identify formulae for the kernels K such that   \begin{eqnarray} c^{(\kappa \kappa ^\prime )}_{st} &=& \int _0^a \alpha _{st}(r)K^{(\alpha ;\kappa \kappa ^\prime )}_{st}(r)+\beta _{st}(r) K_{st}^{(\beta ;\kappa \kappa ^\prime )}(r)\nonumber\\ &&+\,\, \rho _{st}(r)K^{(\rho ;\kappa \kappa ^\prime )}_{st}(r)\,\mathrm{d}r. \end{eqnarray} (9)If the aspherical model is specified in terms of a different set of parameters, equivalent expressions can be found. Thus, an earth model can be obtained by performing a linear inversion upon measured splitting coefficients. A particular attraction of this procedure is the possibility of selecting multiplets that are presumed to have sensitivity at particular depths, and hence to produce models focused on a given region of the Earth (e.g. Koelemeijer et al.2013). 2.2.2 Splitting functions in the real Earth All the foregoing discussion is self-consistent, subject to the fundamental assumption that the (groups of) multiplets are spectrally isolated. However, at least in the context of forward-modelling, it has been conclusively shown that this approach is insufficiently accurate for modern seismology, and instead high-quality seismic calculations must account for interactions throughout the spectrum. Various studies have shown that allowing for cross-coupling in splitting function analysis improves results (e.g. Resovsky & Ritzwoller 1995, 1998; Irving et al.2008; Deuss et al.2013). However, it is not possible to extend the splitting function approach to ever-wider coupling groups: the linearization described in Appendix A4 requires the coupling band-width to be small. If, instead, we choose to work with the full frequency-dependent coupling matrix M(ω), it is unclear whether a splitting-function–style parametrization can be usefully introduced. Where does this leave splitting-function studies? In effect, splitting functions should be regarded as a derived data type or ‘secondary observables’ (e.g. Cara & Lévêque 1987), intended to convey spectral information in a more manageable form. The measurement procedure should therefore be seen as a data-processing operation, or transformation applied to raw seismic data. The particular form of this transformation may be motivated by an out-dated assumption, but it remains a valid transformation to adopt. As such, published measurements of splitting functions may be taken at face value: they need not be regarded as somehow ‘approximate’. The splitting function is defined to be the result of a specified measurement procedure. However, if we are to adopt the isolated-multiplet assumption in the measurement procedure, eq. (9) provides only an approximate relationship between the measured splitting functions and the underlying earth structure. Alternatively, one may take the view that eq. (9) is the defining property of the splitting functions: in this case, the measurement procedure is deficient, since it does not properly account for physical effects in wave propagation. Regardless of the viewpoint adopted, ‘measured’ and ‘predicted’ splitting functions differ, due to approximations within the measurement procedure. Thus, the second key goal of this paper is to assess the impact of this approximation: can models produced using eq. (9) nevertheless be interpreted usefully, or are they too severely contaminated by systematic errors? If fully accurate earth models are to be produced from splitting function data, it is necessary to replace eq. (9) and properly account for multiplet interactions throughout the spectrum. Unfortunately, it is not clear that this can be achieved by simple modifications of the modal depth functions K(r). It seems likely that the most fruitful course would involve obtaining sensitivity kernels via adjoint techniques within a full-coupling framework, allowing the full range of physical effects to be accounted for. However, we do not pursue this idea further within the current paper. Indeed, given the capacity to compute exact sensitivity kernels for spectra using adjoint methods, one may question whether there is value in continuing to adopt the two-stage inversion procedure inherent to splitting function studies. 3 TRUNCATION ERRORS IN FULL-COUPLING SEISMOGRAMS First, we consider the problem of computing seismograms using full-coupling. As we have described, a theoretically complete expression for normal mode seismograms requires summation of an infinite set of seismograms. To make the computations tractable, we must truncate this modal expansion (eq. 1), and only include those modes with eigenfrequencies below some cut-off, ωc. This truncation will inevitably introduce errors into the computed seismograms. How do we choose ωc to ensure that synthetic spectra are accurate below some frequency ωmax? 3.1 The scale of acceptable error In order to address this question, we must somehow define the magnitude of error that we consider acceptable. One way to go about this is to consider the noise in the observed seismograms that we ultimately wish to fit. To allow meaningful analysis, the truncation error should be well below that noise. However, noise levels may vary considerably between different locations and experiments, limiting the generality of this definition. In addition, the earthquake source itself is also not perfectly well known, although it is apparent from inspection of eq. (2) that this uncertainty can be absorbed into the observational data uncertainty, by adding the corresponding variances. Alternatively, given that we typically wish to fit seismograms to infer Earth structure, another definition of an ‘acceptable’ truncation error is that the error should be much smaller than the signal we wish to infer. This is perhaps more straightforward, since it does not depend upon the vagaries of seismic noise. Of course, the signal itself must be above the noise level if we are to make meaningful inferences. We discuss both cases below and show that they can lead to different conclusions depending on the earth model. To quantify these ideas, we will only consider synthetic data calculated for various earth models. Of the three major seismic parameters (S-wave speed, P-wave speed & density), density is known to have the least influence on seismograms, and is a parameter of considerable importance for understanding global geodynamics. We therefore focus on the density signal when assessing whether truncation errors are ‘negligible’. 3.2 Observations using S20RTS We solved eq. (2) exactly, taking all 3-D effects of varying elastic constants and density into account, using the method described in Al-Attar et al. (2012), and computed the same set of seismograms several times varying ωc. In all cases, we used the Global CMT Catalog source mechanism for the 1994 Bolivian event (event code: 060994A) and a global distribution of 129 stations based on those within the Global Seismic Network and IRIS/IDA Seismic Network (see Fig. 1 ). As basis functions, we used eigenfunctions calculated in in the spherically symmetric model PREM (Dziewonski & Anderson 1981). The coupling matrix, explicitly defined in eq. (A10a), is constructed for the shear wave speed model S20RTS (Ritsema et al.1999), which provides values for δ ln Vs with respect to PREM, and where we added scaled compressional wave speed perturbations using δ ln Vp = 0.5 δ ln Vs and scaled density perturbations using δ ln ρ = 0.3 δ ln Vs. These scaling relations are commonly used and are obtained from mineral physics considerations (e.g. Karato 1993). On top of the mantle model, we implemented the simple crustal model from Woodhouse & Dziewonski (1984) and included rotation and ellipticity. In the following, we will consider vertical component spectra. Figure 1. View largeDownload slide Position of the stations (red stars) and the Bolivian event (beach ball) used in this study. Station GRFO is highlighted and shown by a blue star. Figure 1. View largeDownload slide Position of the stations (red stars) and the Bolivian event (beach ball) used in this study. Station GRFO is highlighted and shown by a blue star. We chose ωmax = 3 mHz and used ωc = 3,  4,  5,  and 6 mHz, increasing the number of coupled singlets (toroidals and spheroidals combined) from approximately 2000 to 13 000. We treated the results from ωc = 6 mHz as a set of ‘reference’ synthetic spectra, $$\mathbf {s}_\mathrm{ref}$$, and looked at amplitude spectra of the differences $$\mathbf {s}_{(\omega _c=3\,\mathrm{mHz})}-\mathbf {s}_\mathrm{ref}$$, etc., to provide an indication of the convergence of the full-coupling calculations (Fig. 2). These difference spectra represent the truncation error, and are referred to as such in the figure. It is clear that a truncation at ωc = ωmax can still have substantial errors, but they decrease rapidly as ωc is increased. As far as we know, most published studies that adopt a ‘full-coupling’ strategy have used ωc ≈ ωmax. Figure 2. View largeDownload slide Amplitude spectrum of vertical acceleration at station GRFO (black line) calculated by coupling all modes having frequencies less than ωc = 6 mHz. We refer to this as the ‘reference’ spectrum, sref. Also shown are amplitude spectra of the difference between this reference spectrum and those produced using lower values of ωc (blue lines), which we describe as the ‘truncation error’ for that cut-off. The red line shows the effect of removing any contribution from lateral density variations within the reference spectrum. Note that this is of broadly similar scale to the truncation error at ωc = 3 mHz. The green line shows the effect of perturbing the seismic source with three standard deviations on the reported source uncertainties. Figure 2. View largeDownload slide Amplitude spectrum of vertical acceleration at station GRFO (black line) calculated by coupling all modes having frequencies less than ωc = 6 mHz. We refer to this as the ‘reference’ spectrum, sref. Also shown are amplitude spectra of the difference between this reference spectrum and those produced using lower values of ωc (blue lines), which we describe as the ‘truncation error’ for that cut-off. The red line shows the effect of removing any contribution from lateral density variations within the reference spectrum. Note that this is of broadly similar scale to the truncation error at ωc = 3 mHz. The green line shows the effect of perturbing the seismic source with three standard deviations on the reported source uncertainties. To provide a sense of scale for this truncation error, we repeated the calculation using ωc = 6 mHz, but using a version of S20RTS where all lateral density heterogeneity has been omitted (i.e. with density structure as given by PREM). The amplitude spectrum of $$\mathbf {s}_{\bar{\rho }}-\mathbf {s}_\mathrm{ref}$$ is also shown in Fig. 2 and gives some indication of the magnitude of seismic signal that might be attributable to density structure. Perhaps surprisingly, we observe that this is similar to the magnitude of errors arising from truncation at ωc = ωmax. To estimate the effect of source uncertainties, we make use of those reported in the Global CMT Catalogue, although we note that these are likely to be an under-estimate (Valentine & Trampert 2012). We therefore implement a large source variation by perturbing each parameter by three times the corresponding reported standard deviation, with positive or negative perturbations chosen randomly, and again computed spectra. Given the other effects seen in Fig. 2, the source effect is very modest. We will therefore neglect its contribution in the discussion below. Rather than plotting similar figures for every station around the globe, we sought a method to assess the truncation errors more systematically. We defined a measure of the relative spectral misfit within any given frequency range (ω1, ω2)   \begin{equation} \xi (\mathbf {s},\omega _1,\omega _2) = \frac{\int _{\omega _1}^{\omega _2}\left|\mathbf {s}(\omega )-\mathbf {s}_\mathrm{ref}(\omega )\right|^2\,\mathrm{d}\omega }{\int _{\omega _1}^{\omega _2} \left|\mathbf {s}_\mathrm{ref}(\omega )\right|^2\,\mathrm{d}\omega }. \end{equation} (10)We then plot histograms of ξ evaluated for every station in our data sets (Fig. 3). Again, we compare the reference spectra to those truncated at ωc = 3, 4, and 5 mHz, and to those obtained by using ωc = 6 mHz, but omitting lateral density. If we wish to ensure that the truncation error is consistently smaller than the density signal, we require the corresponding histograms to not overlap. We clearly see that a truncation at 3 mHz will generally lead to errors of a similar scale to the signal anticipated from density. Even for a truncation at 4 mHz some overlap of the distributions occurs for frequencies close to ωmax. It therefore seems necessary to ensure that the truncation frequency is significantly above the maximum frequency of interest (i.e. ωc ≫ ωmax). These results are broadly in accordance with what we might expect based on the knowledge of the coupling-matrix elements. Generally, modes close together in frequency interact most strongly, so truncation errors ought to grow as we approach the cut-off frequency. Figure 3. View largeDownload slide Histograms showing relative spectral misfit, ξ (eq. 10), across all stations. Again, we compare reference spectra (ωc = 6 mHz) to those truncated at lower frequencies (ωc = 3, 4, 5 mHz), with all calculations being performed in the model S20RTS (where we incorporated Vp and density heterogeneity as described in the main text). Results from comparing the same reference spectra to those calculated in a version of S20RTS without any lateral density heterogeneity are also shown. Figure 3. View largeDownload slide Histograms showing relative spectral misfit, ξ (eq. 10), across all stations. Again, we compare reference spectra (ωc = 6 mHz) to those truncated at lower frequencies (ωc = 3, 4, 5 mHz), with all calculations being performed in the model S20RTS (where we incorporated Vp and density heterogeneity as described in the main text). Results from comparing the same reference spectra to those calculated in a version of S20RTS without any lateral density heterogeneity are also shown. To quantify the truncation effects further, we defined a cumulative misfit for the entire synthetic data set   \begin{equation} \psi (\omega ) = \frac{\sum _ i \int _0^\omega |\mathbf {s}^{(i)}(\omega ^\prime )-\mathbf {s}_\mathrm{ref}^{(i)}(\omega ^\prime )|^2\,\mathrm{d}\omega ^\prime }{\sum _i \int _0^\omega |\mathbf {s}_\mathrm{ref}^{(i)}(\omega ^\prime )|^2\,\mathrm{d}\omega ^\prime } \end{equation} (11)where the summation is over all the stations in the data set. We then plot this quantity as a function of frequency for the various cut-off frequencies (Fig. 4), and for the data set where density heterogeneity has been omitted. There are several ways of reading this figure. First, focusing on the curves for the various choices of ωc: suppose we want a precision of (say) 0.1 per cent of the total energy in the signal (this choice should be guided by knowledge of the uncertainties in observed seismic spectra and the seismic source). In this case, truncation at ωc = 3 mHz will produces seismograms accurate up to ωmax = 2 mHz (as shown by a horizontal line on the figure). For a truncation at ωc = 4 mHz, seismograms are accurate up to 3 mHz within the same precision. It appears that we roughly need coupling of 1 mHz higher than the highest frequency of interest in the seismogram. Alternatively, if we are interested in imaging Earth's density structure, we should require truncation errors to be at least an order of magnitude smaller than the density signal in the spectra, and we see roughly that the cut-off frequency should be at 4 mHz for any frequency ωmax ≤ 3 mHz. These values are based on an assumption that Earth's heterogeneity is similar to that present in our implementation of S20RTS. The reader may choose her/his own desired precision and adapt the conclusions accordingly, by direct reading of Fig. 4. Figure 4. View largeDownload slide Cumulative misfits defined by eq. (11), relative to reference spectra computed in S20RTS, for three different cut-off frequencies and for the omission of lateral density heterogeneity. The horizontal grey line denotes a misfit of 0.1  per  cent, as discussed in the main text. Figure 4. View largeDownload slide Cumulative misfits defined by eq. (11), relative to reference spectra computed in S20RTS, for three different cut-off frequencies and for the omission of lateral density heterogeneity. The horizontal grey line denotes a misfit of 0.1  per  cent, as discussed in the main text. 3.3 Model-dependence of the truncation error We expect results to vary depending on the model, and in particular on the power spectrum of heterogeneity. Coupling effects are somewhat analogous to scattering, and it is well-known that short-wavelength heterogeneity promotes multiple-scattering effects. We might reasonably expect to observe similar effects in coupling calculations, with stronger short-wavelength structure necessitating relatively higher values of ωc. S20RTS has modest amplitudes and a ‘red’ spectrum, dominated by long-wavelength structure. The model resolution is far below what the parametrization of the model would allow, ranging from 2000–4000 km horizontally to 250–750 km vertically (Ritsema et al.2004). It is therefore important to investigate how this is dictating the findings in the previous subsection. To provide some insight into this question, we generated a random model using the same parametrization as S20RTS—a spherical harmonic expansion to degree 20 laterally and 21 splines vertically. For δ ln Vs we drew random numbers uniformly in the range [−0.01, 0.01] giving an rms-amplitude of about 12 per cent at all depths. We continue to assume that Vp and ρ are scaled versions of Vs, as with S20RTS. The rms-amplitude of this random model is between 5 and 38 times stronger than S20RTS depending on depth. This produces a model of true degree-20 lateral resolution, a white spectrum, and much higher vertical resolution. We repeated all the full-coupling calculations already described, using this model. Again, we produce histograms of relative spectral misfit (Fig. 5) and cumulative misfit (Fig. 6 ). From Fig. 5, we see that now the truncation effects are much stronger and overall we need to couple more modes than was necessary in S20RTS. In other words, more power in the structural model, and especially at short-wavelengths, necessitates higher values of ωc. It is also obvious that the contribution of density to the seismograms is much lower in the random model than in S20RTS, despite the fact that the rms-amplitude of the model is much larger. This is because density imprints itself differently into the seismograms, compared to the model with a red spectrum. Figure 5. View largeDownload slide Histograms of relative spectral misfit across all stations defined by eq. (10) for a random model. Figure 5. View largeDownload slide Histograms of relative spectral misfit across all stations defined by eq. (10) for a random model. Figure 6. View largeDownload slide Cumulative misfits defined by eq. (11) for a randomly generated, spectrally white earth model. Compare Fig. 4. Figure 6. View largeDownload slide Cumulative misfits defined by eq. (11) for a randomly generated, spectrally white earth model. Compare Fig. 4. Clearly the relation between ωc, ωmax and the earth model is not straightforward. Looking at cumulative misfits (Fig. 6), we see that white models appear to yield larger truncation errors than their red counterparts, as expected. Again, considering the data sets that include density heterogeneity, we can imagine seeking a precision of 0.1 per cent of the total energy in the signal. A truncation at ωc = 3 mHz now only produces seismograms accurate up to 1 mHz. For a truncation at ωc = 4 mHz, seismograms are accurate up to around 2 mHz within the same precision. It appears that we now need coupling to be at least 2 mHz higher than the highest frequency of interest in the seismograms. If, however, we define the acceptable truncation error as being an order of magnitude smaller than the signal arising from density heterogeneity, we now see that coupling up to 5 mHz is required for ωmax ≤ 3 mHz. Again it is straightforward to adapt the conclusions for any other desired precision. 3.4 Summary The key conclusions we draw from these experiments may be summarized: Accurate full-coupling spectra require the cut-off frequency ωc to be significantly above the maximum frequency present in the seismograms, ωmax. To our knowledge, this criterion has not been met in previously published studies based on full-coupling. Accurate seismograms up to 3 mHz in earth-like models (having a red spectrum) require truncation at no less than 4 mHz. Generally truncation errors are stronger at higher frequencies, as might be expected. Effects seem more pronounced as models gain more power, especially at shorter wavelengths. In models with white spectra, accurate seismograms up to 3 mHz require a truncation level of at least 5 mHz. Overall, these observations are all broadly consistent with expectations based on theoretical arguments. However, the effects are now quantified, and we have obtained rules-of-thumb that can be applied to ensure acceptable accuracy when performing calculations. 4 THE INTERPRETABILITY OF SPLITTING FUNCTION MEASUREMENTS We now turn to consider splitting functions, and the extent to which they can provide useful constraints upon Earth structure. Catalogues of splitting function measurements (e.g. Giardini et al.1987; He & Tromp 1996; Resovsky & Ritzwoller 1998; Durek & Romanowicz. 1999; Masters et al.2000; Deuss et al.2011, 2013) still represent the main ‘derived’ data set for long-period seismology. Measured centre frequencies (i.e. ω0 in eq. 5) are taken as representative of the Earth's spherical average structure (e.g. de Wit et al.2014) and splitting function coefficients feature as a constraint in many current global tomography models (e.g. Ritsema et al.2011; Moulik & Ekström 2016). Splitting functions have also been used for targeted studies investigating specific aspects of Earth structure, since one can examine the kernels from eq. (9) and identify modes which have the desired sensitivity (see Fig. 7). This has particularly been exploited to investigate deep-earth structure, including inner core studies (beginning with Woodhouse et al.1986), targeted lower mantle studies (e.g. Koelemeijer et al.2013, 2017) and density inferences (e.g. Ishii & Tromp 1999; Resovsky & Trampert 2003). Figure 7. View largeDownload slide Self-coupling sensitivity kernels $$K_{00}^{(\kappa \kappa )}$$ (as in eq. 9) corresponding to the angular order s = 0 splitting coefficient for our selected modes. The solid red line denotes Vp sensitivity, the solid black line Vs sensitivity and the dotted black line gives density sensitivity. The modes are ordered somewhat arbitrarily, with a general trend from mostly upper-mantle–sensitive to mostly lower-mantle–sensitive. Figure 7. View largeDownload slide Self-coupling sensitivity kernels $$K_{00}^{(\kappa \kappa )}$$ (as in eq. 9) corresponding to the angular order s = 0 splitting coefficient for our selected modes. The solid red line denotes Vp sensitivity, the solid black line Vs sensitivity and the dotted black line gives density sensitivity. The modes are ordered somewhat arbitrarily, with a general trend from mostly upper-mantle–sensitive to mostly lower-mantle–sensitive. It is common for splitting function catalogues to publish uncertainty estimates on the recovered splitting coefficients and centre-frequencies. In principle, these can be obtained directly from the Bayesian inversion framework (e.g. Tarantola & Valette 1982). However, it is not clear how to define a meaningful prior for the splitting function inversion, and in any case such analysis would neglect the inherent non-linearity of the measurement process. This point has previously been made by Resovsky & Ritzwoller (1998). After a careful analysis of all potential contributions to the uncertainty, they estimated combined uncertainties using regressions performed on synthetic data with simulated theoretical errors and data uncertainties. More recently, it has been common to use some form of boot-strapping procedure to investigate the degree of constraint upon results (e.g. Deuss et al.2013). This uncertainty estimate therefore reflects any apparent inconsistency between different parts of the data set: it may not reflect any systematic effects inherent to the measurement procedure or its implementation. Various studies have highlighted that considering interactions between adjacent modes can alter splitting-function results significantly (e.g. Resovsky & Ritzwoller 1995; Irving et al.2008) and group-coupling, where appropriate, is now the norm for various frequency bands (e.g. Deuss et al.2013). Studies have also investigated the extent to which splitting function inversions are well-constrained, typically focusing on aspects of the linearized inversion scheme, such as regularization, parametrization and starting model (e.g. Resovsky & Ritzwoller 1999; Kuo & Romanowicz 2002; Pachhai et al.2016). However, we are not aware of any previous studies that have taken a more holistic approach, to investigate whether splitting functions provide a meaningful representation of Earth structure that can be interpreted by eq. (9). As already discussed, a potential issue arises because the measurement procedure is inherently based upon isolated-group approximations (either self- or group-coupling), although their interpretation (via eq. 9) is not. Our results in Section 3, and those of earlier work (e.g. Deuss & Woodhouse 2001; Al-Attar et al.2012; Yang & Tromp 2015), all indicate that isolated-group seismograms provide a poor approximation to the true physics of wave propagation in the earth. It therefore seems distinctly possible that structural models estimated via isolated-group synthetics may inherit these inaccuracies. We address this by performing a straightforward synthetic test. For any earth model, we compute ‘predicted’ self-coupling splitting function coefficients using eq. (9). We also compute accurate synthetic seismograms, using full-coupling theory and the results of the previous section. From these seismograms, we then obtain ‘measured’ splitting functions $$c^{(\kappa \kappa )}_{st}$$, following the same procedure as is used for real data. We then assess the agreement between measured and predicted splitting function values. Ideally, they should agree to within the stated uncertainty for the measurement procedure. If this is the case, we can conclude that the measured splitting functions are interpretable using eq. (9), and the measurement uncertainties can be propagated through this expression to provide uncertainties on structural parameters. If measured and predicted splitting functions differ by more than the measurement uncertainty—but show no evidence of any systematic bias—it may be reasonable to continue to apply eq. (9) for interpretation using increased uncertainty levels to reflect the theory error (as per Tarantola & Valette (1982)). However, if systematic biases are found, no meaningful interpretation can be made. 4.1 Measurement of synthetic splitting functions As stated in the detailed analysis of Resovsky & Ritzwoller (1998), splitting function measurements depend on theoretical errors, regularization, noise in the data and event-station distributions. The final splitting function depends non-linearly on all of these. We want to quantify as precisely as possible the theoretical errors, and so we need to specify the other contributions as accurately as possible. To do this, all details matter. We therefore generated a synthetic full-coupling data set and measured splitting functions following the approach of Deuss et al. (2013). We considered these new spectra to be the observed data, referred to as above as $$\mathbf {s}_\mathrm{ref}$$. We will only focus on the self-coupling splitting function η(κκ) (eq. 7), although measurements in certain parts of the spectrum employed group-coupling. We used the same event and station distributions and method (forward theory, cut-off frequencies and regularization) as those in Deuss et al. (2013). The measured splitting coefficients $$c^{(\kappa \kappa )}_{st}$$ can then directly be compared to the expected coefficients for model S20RTS (as obtained from eq. 9). Two factors contribute to any difference between measured and expected coefficients. First, differences may be the result of the approximations in the measurement procedure: this is the effect we wish to address. However, they may also arise as a result of the regularization applied within the non-linear optimization procedure. To separate the two contributions, we also calculated S20RTS seismograms strictly using the self-coupling approximation, and again measured splitting functions from these. In this case, any difference between measured and observed coefficients can be attributed solely to regularization, as data and measurement procedure now contain the same physics. Assuming that all contributions are Gaussian distributed, the difference between the covariances corresponding to the two sets of splitting function then isolates the theoretical error. We focused on a selection of modes, which were also measured by Deuss et al. (2013), namely: 0S2, 0S5, 0S6, 0S7, 1S4, 1S7, 1S8, 1S9, 1S10, 2S6, 2S12, 2S13, 5S3. Fig. 7 shows some of the kernels K, which we used in eq. (9) to predict the splitting functions we ought to retrieve if the measurements were perfect. Some of the chosen modes have sensitivity mainly in the lower mantle, some mainly in the upper mantle, and some are sensitive throughout the mantle. Our discussion will concentrate on coefficients of angular order 2 and 4, which carry most of the splitting function energy for the modes we selected to analyse. This will also avoid the problem of neglected higher degree structure as the measurements are done to higher orders (Resovsky & Ritzwoller 1998). All our synthetic experiments are noise-free because we want to isolate the theoretical error without a particular noise model (which would be difficult to choose objectively) propagating non-linearly into the results via the regularization. Assuming that splitting function measurements are secondary data to be taken at face value, we will put our synthetic measurements into context by using actual uncertainties from Deuss et al. (2013) obtained by boot-strapping real earthquake data and thought to represent only the data noise present in the Deuss catalogue. 4.2 Observations from synthetic splitting functions The coefficients of the splitting function measurements for mode 0S6, a lower mantle mode, are depicted in Fig. 8. When self-coupling seismograms are used, the measured coefficients are indistinguishable from those calculated for S20RTS via the kernels, suggesting that the regularization used for this mode is not affecting the measurement. We therefore only show the calculated coefficients in the figure, but not those measured from self-coupling seismograms. There is however a difference with those measured using the full-coupling spectra as data. To calibrate this difference against the density signal, we also computed predicted splitting coefficients for a version of S20RTS without lateral density heterogeneity. While there are differences between the three sets of splitting coefficients, these are small compared to the magnitude of the coefficients. It is therefore instructive to plot differences. Figure 8. View largeDownload slide Observed splitting coefficients for mode 0S6 for s = 2 and s = 4 (top panels). The full black circles represent measurements made on synthetic spectra obtained using full-coupling. The red open circles are predicted splitting coefficients calculated using eq. (9), and correspond to those you would expect if the measurements were perfect. The grey diamonds represent splitting coefficients calculated for mantle model S20RTS without lateral density variations. Bottom panels: difference between those measured from full-coupling spectra and predicted coefficients (black circles) and difference between predicted coefficients without and with density (grey diamonds). The vertical bars correspond to error bars inferred by Deuss et al. (2013) for this mode. Figure 8. View largeDownload slide Observed splitting coefficients for mode 0S6 for s = 2 and s = 4 (top panels). The full black circles represent measurements made on synthetic spectra obtained using full-coupling. The red open circles are predicted splitting coefficients calculated using eq. (9), and correspond to those you would expect if the measurements were perfect. The grey diamonds represent splitting coefficients calculated for mantle model S20RTS without lateral density variations. Bottom panels: difference between those measured from full-coupling spectra and predicted coefficients (black circles) and difference between predicted coefficients without and with density (grey diamonds). The vertical bars correspond to error bars inferred by Deuss et al. (2013) for this mode. In the bottom panels of Fig. 8 we plot the difference between the splitting functions predicted for S20RTS, and those measured from full-coupling synthetics (black line). Each coefficient is plotted with error-bars, representing the uncertainty reported for this mode by Deuss et al. (2013). We see that most coefficients touch the zero-line when the error-bars are taken into account (certainly when 2 standard deviations or the 95 per cent confidence level is considered), implying that predicted and measured splitting functions agree to within the measurement uncertainty. Thus, for this mode, any systematic errors that arise from the use of self-coupling are negligible compared to the effect of noise within the data. We also plot the difference between splitting functions predicted in S20RTS, and our version without density (grey line). Again, we plot error-bars representing the uncertainty reported by Deuss et al. (2013), and find that many of these also touch the zero-line and certainly within 2 standard deviations. This suggests that the signal from density lies at, or below the noise level of measured splitting functions. As a result, it is not straightforward to extract meaningful information about density from splitting functions for this mode. In addition, we observe that the signal attributable to lateral density heterogeneity appears to be similar in magnitude to the difference between observed and predicted splitting coefficients. Thus, even if noise levels could be reduced significantly, it would appear that current splitting function measurement procedures are insufficiently accurate to make use of the density signal. Similar results are seen in Fig. 9, which shows results for the upper mantle mode 2S12. In this example, we see differences which are more significant relative to uncertainties. Errors due to the use of self-coupling in the measurement procedure are now markedly greater than the quoted uncertainty, which accounts only for effects due to noise in the data, even within 2 standard deviations for some coefficients. This is likely to affect any interpretation based on eq. 9. Although lateral density heterogeneity has a signal that is now also substantially greater than the stated noise level, it remains comparable in size to the theoretical error in splitting function observations. Again, it seems unlikely that meaningful density inferences can be made from this mode using standard methods. These two examples suggest that both density signal, and measurement error, become more significant for modes with shallow sensitivity. To see if this holds in general, we must examine more modes. Figure 9. View largeDownload slide As Fig. 8, but for mode 2S12. Figure 9. View largeDownload slide As Fig. 8, but for mode 2S12. Comparing the full suite of splitting function coefficients measured from full-coupling synthetics to those predicted using eq. (9), we find a very strong correlation (with an R-value of 0.99), with no obvious evidence of any systematic trends. Splitting functions depend non-linearly on the spectra, and a more sophisticated error analysis is required before the possibility of any systematic effects can be fully eliminated. However, our results suggest that a reasonable starting point may be to assume that theoretical measurement errors can be modelled by a zero-mean Gaussian distribution, allowing them to be straightforwardly incorporated into measurement uncertainties. To interpret our observations, it is then most instructive to show the reduced χ2 misfits for all measured coefficients of a given splitting function with respect to the coefficients predicted for S20RTS using eq. (9)   \begin{equation} \chi ^2 = \frac{1}{N} \sum _{s,t} \frac{\left(c^\mathrm{meas.}_{st}-c^\mathrm{pred.}_{st}\right)^2}{\sigma ^2_{st}} \end{equation} (12)where σst are the uncertainties from Deuss et al. (2013) for the given mode and N represents the number of coefficients in the corresponding splitting function. This statistical measure allows results to be summarized for easy interpretation; χ2 = 1 indicates that the difference between measured and predicted values is, on average, of similar size to the uncertainties. First, we consider the splitting function measurements made upon synthetic seismograms computed using full-coupling (see Fig. 10 ). For most modes, the value of χ2 is well above one, indicating that observations differ from predicted values by an amount significantly greater than currently published data uncertainties for splitting coefficients. In other words, the Deuss catalogue underestimates the true uncertainty on splitting coefficients. This in turn may impact upon the confidence that can be placed in models derived partly or fully from these splitting function data. Figure 10. View largeDownload slide Reduced χ2 misfit defined by eq. (12) for measurements on full- (black squares) and self-coupling (open squares) seismograms using all parameters as determined by Deuss et al. (2013). Figure 10. View largeDownload slide Reduced χ2 misfit defined by eq. (12) for measurements on full- (black squares) and self-coupling (open squares) seismograms using all parameters as determined by Deuss et al. (2013). As has already been discussed, discrepancies between measured and predicted splitting coefficients may arise due to reliance upon approximations in the measurement procedure, or as a consequence of regularization. To compare these effects, Fig. 10 also shows χ2 for splitting functions measured from synthetic seismograms generated using self-coupling. In this case, seismogram generation and splitting function measurement employ a consistent theoretical framework. Thus, any differences here arise solely from regularization. For several modes, we see that χ2 obtained from self-coupling synthetics is similar to that for full-coupling synthetics. This suggests that the discrepancy between prediction and observation in these cases is mostly attributable to regularization. It is possible that reducing the damping applied during the measurement procedure would improve the agreement between predictions and observations. However, this might also be expected to increase the sensitivity to data noise, resulting in increased error-bars on the recovered coefficients (Resovsky & Ritzwoller 1998). It is therefore not obvious that a net benefit would be gained; re-analysis of real data to investigate this rigorously is beyond the scope of the present study. Assuming that all contributions to the uncertainties are Gaussian distributed, covariances will simply add up, in other words the total posterior covariance in actually measured coefficients may be expressed as $$\mathbf {C}_\mathrm{tot}=\mathbf {C}_\mathrm{th}+\mathbf {C}_\mathrm{reg}+\mathbf {C}_\mathrm{d}$$, where $$\mathbf {C}_\mathrm{th}$$ is the theoretical uncertainty we wish to quantify, $$\mathbf {C}_\mathrm{reg}$$ is due to regularization in the regression and $$\mathbf {C}_\mathrm{d}=\sigma ^2_d \mathbf {I}$$ corresponds to data uncertainty including the influence of the event-station distribution. This decomposition is similar to that used in Resovsky & Ritzwoller (1998). The uncertainties quoted in Deuss et al. (2013) were determined using boot-strapping: randomly eliminating stations, and also whole events. This was done for a fixed regularization. Thus, these uncertainties therefore represent the data covariance, $$\mathbf {C}_\mathrm{d}$$, including the particular event-station distribution used to create the Deuss catalogue. Our measurements based on self-coupling seismograms (see Fig. 10) do not reflect the calculated splitting functions perfectly. Because these seismograms are noise free, the discrepancies may be a measure of the regularization noise, which can be expressed as a multiple of the data noise, $$\mathbf {C}_\mathrm{reg}=\alpha \sigma ^2_d \mathbf {I}$$. Values of α can be directly read from Fig. 10. In the same figure we also display the result for measurements on full-coupling seismograms. These represent the discrepancies due to regularization and the theoretical error combined, which correspond then to $$\mathbf {C}_\mathrm{th}+\mathbf {C}_\mathrm{reg}=\beta \sigma ^2_d \mathbf {I}$$. Hence β − α is a measure of $$\mathbf {C}_\mathrm{th}$$ expressed in units of the covariances in the Deuss catalogue. We see that $$\mathbf {C}_\mathrm{th}$$ and $$\mathbf {C}_\mathrm{reg}$$ changes for every mode, but as a quantitative indication: for the sample of modes we considered, on average β − α = 1.3 and α = 2.0. Thus for the Deuss catalogue, the theoretical measurement error, is on average 1.3 times the quoted variance and the regularization error is on average 2 times the quoted variance, resulting in $$\mathbf {C}_\mathrm{tot}=(1+1.3+2.0)\sigma ^2_d\mathbf {I}=4.3\sigma ^2_d\mathbf {I}$$. We therefore recommend users should multiply the published uncertainty in the Deuss catalogue by a factor of 2 to use a more representative number for the true measurement uncertainty. While the quoted numbers reflect the contribution to uncertainty in the Deuss catalogue, can they be used for other catalogues? The answer is no, because regularization errors and the data uncertainty translate differently into coefficient uncertainty, depending intimately on all the details on the measurement procedure (event-station distribution, regularization, number of iterations, cut-off frequencies, …). Because we do not have this information available for other catalogues, we cannot give detailed results for those. We can however say Resovsky & Ritzwoller (1998) made a real effort to quantify the true uncertainties given the computational restrictions of 30 years ago. To put these results into perspective, it is instructive to look at potential information on the Earth's structure contained in the splitting functions. We have already described how we generate a version of S20RTS omitting any lateral density heterogeneity, and use this to estimate the magnitude of signal attributable to density. Similarly, we can generate seismograms in versions of S20RTS that assume there is no lateral Vs heterogeneity (having only lateral variations in Vp and ρ), or no lateral Vp structure (only Vs and ρ). In each case, we can compute expected splitting coefficients using eq. (9), and obtain χ2 values as in eq. (12). When we plot these in Fig. 11, we see that the signal of Vs in the splitting functions is 2-3 orders of magnitude larger than the observational uncertainties for most modes, while that of Vp is on average one order of magnitude larger and ρ mostly of the same order as the uncertainties. This is interesting in itself, as it shows that most modes do not contain any significant information on density, regardless of whether we can accurately measure their splitting function or not. The reason for this might be that wave speeds mostly affect the phase of a seismic spectrum, whereas density mostly affects its amplitude. When we further consider our findings that the measured data variances should be multiplied by a factor of 4.3 to represent the total measurement errors, we find that on average the reduced χ2 will be around 100 for δ ln Vs, around 4 for δ ln Vp and around 1 for δ ln ρ. Figure 11. View largeDownload slide Reduced χ2 misfit defined by eq. (12) for model S20RTS, except with δ ln Vs = 0 (black squares), δ ln Vp = 0 (open squares) and δ ln ρ = 0 (grey diamonds). This provides an indication of the average strength of the respective fields within the splitting functions for the 13 modes analysed. Figure 11. View largeDownload slide Reduced χ2 misfit defined by eq. (12) for model S20RTS, except with δ ln Vs = 0 (black squares), δ ln Vp = 0 (open squares) and δ ln ρ = 0 (grey diamonds). This provides an indication of the average strength of the respective fields within the splitting functions for the 13 modes analysed. Most studies infer mantle structure using these splitting function in a least-squares inversion (e.g. Ritsema et al.1999; Ishii & Tromp 1999; Ritsema et al.2011; Moulik & Ekström 2016; Koelemeijer et al.2017). Given that the Vs signal is well above the uncertainty of the splitting functions, shear wave mantle models using such data should be fairly robust. However, because the density signal is of the level of the measurement uncertainty, its robustness was justifiably questioned in the past (e.g. Resovsky & Ritzwoller 1999; Kuo & Romanowicz 2002). Only by moving away from standard techniques and using full sampling, can density be inferred with a meaningful uncertainty (Resovsky & Trampert 2003). The latter study showed that the density model had on average an uncertainty of about 50 per cent of its magnitude, or 100 per cent within a 95 per cent confidence level. This generated lively debates, and for instance, the question of buoyant versus heavy large low-shear wave provinces is still waiting to be settled, although the heavy side (Ishii & Tromp 1999; Resovsky & Trampert 2003; Trampert et al.2004; Mosca et al.2012; Moulik & Ekström 2016) got a recent boost with an independent study inverting tidal data (Lau et al.2017). Given that the Deuss catalogue uses far more seismograms to infer splitting functions than any previous catalogue, and we now have a good representations of the true measurement uncertainty, constructing a mantle model using a sampling algorithm is well worth repeating. Alternatively, direct spectra inversions using full coupling theory, where density has a higher signal-to-noise ratio, are becoming computationally feasible. 4.3 Summary Again, our results can be summarized as: We confirm that there is a significant theoretical measurement error in currently published splitting functions. There is no obvious bias, but this theoretical error can, to a first approximation, be assumed to be random noise. There is also a significant contribution from regularization to the uncertainty in splitting functions. Taking the total uncertainty into account, currently published splitting functions contain information on Vs with a signal-to-noise ratio of about 100, information on Vp with a signal-to-noise ratio of about 4, and information on ρ with a signal-to-noise ratio of about 1. Since the density signal is small in currently published splitting functions, advanced techniques should be employed to infer density structure. The signal from Vs is very strong, and its inference is unlikely to be significantly affected by the increased measurement uncertainty in currently published splitting functions. Thus, we do not believe our results should affect published models of Vs structure that may be derived these splitting function measurements. These results are specific to the Deuss catalogue (Deuss et al.2013). Using other catalogues requires repeating the work in this study if authors wish to infer density, the exception possibly being the Colorado catalogue (Resovsky & Ritzwoller 1998). We also remark that the density signal may be more effectively isolated by working directly upon seismic spectra, where the observational uncertainty is expected to be smaller: the results of Section 3 suggest that this should be feasible. 5 CONCLUSIONS Previous studies have shown that self- and group-coupling approximations ought to be regarded as obsolete in the context of computing seismograms in 3-D Earth models. By systematically increasing the number of modes in full-coupling calculations, we have demonstrated that accurate full-coupling requires inclusion of modes significantly beyond the highest frequency range of interest in the seismograms. For seismograms to be accurate up to 3 mHz in typical models, coupling to 5 mHz is preferable; 4 mHz may be sufficient for some applications if the Earth has a structure similar to S20RTS. Our experimental evidence appears to corroborate theoretical expectations that truncation errors become larger as short-wavelength structure becomes more dominant. Self-coupling splitting functions, often used in the construction of 3-D Earth models, are themselves inherently based on the self- or group-coupling approximations, and may therefore contain systematic measurement errors. By measuring self-coupling splitting functions on exact full-coupling seismograms, we showed that the recovered splitting function does not match theoretical predictions. There is also a significant regularization contribution to the recovered splitting functions. There is no visible apparent bias from any of these contributions, but recently published (Deuss et al.2013) standard deviations of data uncertainties for splitting coefficients need to be increased by a factor of around 2 to represent the total error within the measurement process. In that case, our results show that the information on the 3-D density structure contained in the measured splitting functions has the same magnitude as the total uncertainty. The information on the variation of compressional-wave speed is about 4 times that of the uncertainty and the information on the shear-wave speed has a signal-to-noise ratio of about 100. To go beyond a statistical treatment of the theoretical error, one route forward may be to simply avoid the use of splitting functions: with modern computational resources, it may be feasible to infer structural information directly from spectra via full-coupling calculations. Alternatively, it may also be possible to continue to use splitting functions, measured using self- or group-coupling approximations, and instead replace the kernels within eq. (9) by those that account for these measuring approximations correctly (e.g. via adjoint theory). However, until such studies are undertaken, we must assume markedly larger uncertainties than reported for measured splitting coefficients. Given this, our modelling suggests that Earth's density signal lies at or below the presumed noise level. In consequence, there is no prospect of recovering robust estimates of density structure. It therefore follows that currently published density models ought to be interpreted with great caution. ACKNOWLEDGEMENTS We are grateful for constructive reviews by Philippe Lognonné and Joseph Resovsky. The research leading to these results has received funding from the European Research Council (ERC) under the European Union's Seventh Framework Programme (FP/2007-2013) grant agreement number 320639 (iGEO) and under the European Union's Horizon 2020 research and innovation programme grant agreement number 681535 (ATUNE). The code used to perform the full-coupling calculations makes extensive use of routines developed by John Woodhouse, and we thank him for allowing us to use them. REFERENCES Al-Attar D., Crawford O., 2016. Particle relabelling transformations in elastodynamics, Geophys. J. Int.  205 575– 593. https://doi.org/10.1093/gji/ggw032 Google Scholar CrossRef Search ADS   Al-Attar D., Woodhouse J., Deuss A., 2012. Calculation of normal mode spectra in laterally heterogeneous earth models using an iterative direct solution method, Geophys. J. Int.  189 1038– 1046. https://doi.org/10.1111/j.1365-246X.2012.05406.x Google Scholar CrossRef Search ADS   Cara M., Lévêque J., 1987. Waveform inversion using secondary observables, Geophys. Res. Lett.  14 1046– 1049. https://doi.org/10.1029/GL014i010p01046 Google Scholar CrossRef Search ADS   Dahlen F., 1968. The normal modes of a rotating, elliptical Earth, Geophys. J. R. astr. Soc.  16 329– 367. https://doi.org/10.1111/j.1365-246X.1968.tb00229.x Google Scholar CrossRef Search ADS   Dahlen F., 1969. The normal modes of a rotating, elliptical Earth—II. Near-resonance multiplet coupling, Geophys. J. R. astr. Soc.  18 397– 436. https://doi.org/10.1111/j.1365-246X.1969.tb03576.x Google Scholar CrossRef Search ADS   Dahlen T., Tromp J., 1998. Theoretical Global Seismology  Princeton University Press. de Wit R., Käufl P., Valentine A., Trampert J., 2014. Bayesian inversion of free oscillations for Earth's radial (an)elastic structure, Phys. Earth planet. Inter.  237 1– 17. https://doi.org/10.1016/j.pepi.2014.09.004 Google Scholar CrossRef Search ADS   Deuss A., Woodhouse J., 2001. Theoretical free-oscillation spectra: the importance of wide band coupling, Geophys. J. Int.  146 833– 842. https://doi.org/10.1046/j.1365-246X.2001.00502.x Google Scholar CrossRef Search ADS   Deuss A., Woodhouse J., 2004. Iteration method to determine the eigenvalues and eigenvectors of a target multiplet including full mode coupling, Geophys. J. Int.  159 326– 332. https://doi.org/10.1111/j.1365-246X.2004.02399.x Google Scholar CrossRef Search ADS   Deuss A., Ritsema J., van Heijst H., 2011. Splitting function measurements for Earth's longest period normal modes using recent large earthquakes, Geophys. Res. Lett.  38 L04303, doi:10.1029/2010GL046115. https://doi.org/10.1029/2010GL046115 Google Scholar CrossRef Search ADS   Deuss A., Ritsema J., van Heijst H., 2013. A new catalogue of normal-mode splitting function measurements up to 10 mHz, Geophys. J. Int.  193 920– 937. https://doi.org/10.1093/gji/ggt010 Google Scholar CrossRef Search ADS   Durek J.J., Romanowicz B., 1999. Inner core anisotropy inferred by direct inversion of normal mode spectra, Geophys. J. Int.  139 599– 622. https://doi.org/10.1046/j.1365-246x.1999.00961.x Google Scholar CrossRef Search ADS   Dziewonski A., Anderson D., 1981. Preliminary reference Earth model, Phys. Earth planet. Inter.  25 297– 356. https://doi.org/10.1016/0031-9201(81)90046-7 Google Scholar CrossRef Search ADS   Gharti H., Tromp J., 2017. A spectral-infinite-element solution of Poisson's equation: an application to self gravity , arXiv:1706.00855v1. Giardini D., Li X.-D., Woodhouse J., 1987. Three-dimensional structure of the Earth from splitting in free-oscillation spectra, Nature  325 405– 411. https://doi.org/10.1038/325405a0 Google Scholar CrossRef Search ADS   Giardini D., Li X.-D., Woodhouse J., 1988. Splitting functions of long-period normal modes of the Earth, J. geophys. Res.  93 13 716– 13 742. https://doi.org/10.1029/JB093iB11p13716 Google Scholar CrossRef Search ADS   Gilbert F., 1970. Excitation of the normal modes of the Earth by earthquake sources, Geophys. J. R. astr. Soc.  22 223– 226. https://doi.org/10.1111/j.1365-246X.1971.tb03593.x Google Scholar CrossRef Search ADS   Hara T., Tsuboi S., Geller R., 1991. Inversion for laterally heterogeneous earth structure using a laterally heterogeneous starting model: preliminary results, Geophys. J. Int.  104 523– 540. https://doi.org/10.1111/j.1365-246X.1991.tb05699.x Google Scholar CrossRef Search ADS   Hara T., Tsuboi S., Geller R., 1993. Inversion for laterally heterogeneous upper mantle S-wave velocity structure using iterative waveform inversion, Geophys. J. Int.  115 667– 698. https://doi.org/10.1111/j.1365-246X.1993.tb01487.x Google Scholar CrossRef Search ADS   He X., Tromp J., 1996. Normal-mode constraints on the structure of the Earth, J. geophys. Res.  101 20 053– 20 082. https://doi.org/10.1029/96JB01783 Google Scholar CrossRef Search ADS   Irving J., Deuss A., Andrews J., 2008. Wide-band coupling of Earth's normal modes due to anisotropic inner core structure, Geophys. J. Int.  174 919– 929. https://doi.org/10.1111/j.1365-246X.2008.03824.x Google Scholar CrossRef Search ADS   Ishii M., Tromp J., 1999. Normal-mode and free-air gravity constraints on lateral variations in velocity and density of Earth's mantle, Science  285 1231– 1236. https://doi.org/10.1126/science.285.5431.1231 Google Scholar CrossRef Search ADS PubMed  Karato S.-i., 1993. Importance of anelasticity in the interpretation of seismic tomography, Geophys. Res. Lett.  15 1623– 1626. https://doi.org/10.1029/93GL01767 Google Scholar CrossRef Search ADS   Koelemeijer P., Deuss A., Ritsema J., 2013. Observations of core-mantle boundary Stoneley modes, Geophys. Res. Lett.  40 2557– 2561. https://doi.org/10.1002/grl.50514 Google Scholar CrossRef Search ADS   Koelemeijer P., Deuss A., Ritsema J., 2017. Density structure of Earth's lowermost mantle from Stoneley mode splitting observations, Nat. Commun.  8 15241. https://doi.org/10.1038/ncomms15241 Google Scholar CrossRef Search ADS PubMed  Komatitsch D., Tromp J., 2002. Spectral-element simulations of global seismic wave propagation—I. Validation, Geophys. J. Int.  149 390– 412. https://doi.org/10.1046/j.1365-246X.2002.01653.x Google Scholar CrossRef Search ADS   Kuo C., Romanowicz B., 2002. On the resolution of density anomalies in the Earth's mantle using spectral fitting of normal-mode data, Geophys. J. Int.  150 162– 179. https://doi.org/10.1046/j.1365-246X.2002.01698.x Google Scholar CrossRef Search ADS   Lau H., Yang H.-Y., Tromp J., Mitrovica J., Latychev K., Al-Attar D., 2015. A normal mode treatment of semi-diurnal body tides on an aspherical, rotating and anelastic Earth, Geophys. J. Int.  202 1392– 1406. https://doi.org/10.1093/gji/ggv227 Google Scholar CrossRef Search ADS   Lau H., Mitrovica J., Davies J., Tromp J., Yang H.-Y., Al-Attar D., 2017. Tidal tomography constrains Earth's deep-mantle buoyancy, Nature  551 321– 326. https://doi.org/10.1038/nature24452 Google Scholar CrossRef Search ADS PubMed  Lay T., 2007. Deep earth structure—lower mantle and D″, in Treatise on Geophysics  Vol. 1, pp. 620– 654, eds Romanowicz B., Dziewoński A., Elsevier. Lognonné P., 1991. Normal modes and seismograms in an anelastic rotating Earth, J. geophys. Res.  96 20 309– 20 319. https://doi.org/10.1029/91JB00420 Google Scholar CrossRef Search ADS   Masters G., Laske G., Gilbert F., 2000. Matrix autoregressive analysis of free-oscillation coupling and splitting, Geophys. J. Int.  143 478– 489. https://doi.org/10.1046/j.1365-246X.2000.01261.x Google Scholar CrossRef Search ADS   Masters G., Barmine M., Kientz S., 2011. Mineos User Manual  Computational Infrastructure for Geodynamics, v. 1.0.2. Mosca I., Cobden L., Deuss A., Ritsema J., Trampert J., 2012. Seismic and mineralogical structures of the lower mantle from probabilistic tomography, J. geophys. Res.  117 B06304, doi:10.1029/2011JB008851. https://doi.org/10.1029/2011JB008851 Google Scholar CrossRef Search ADS   Moulik P., Ekström G., 2016. The relationships between large-scale variations in shear velocity, density, and compressional velocity in the Earth's mantle, J. geophys. Res.  121 2737– 2771. https://doi.org/10.1002/2015JB012679 Google Scholar CrossRef Search ADS   Pachhai S., Tkalčić H., Masters G., 2016. Estimation of splitting functions from Earth's normal mode spectra using the neighbourhood algorithm, Geophys. J. Int.  204 111– 126. https://doi.org/10.1093/gji/ggv414 Google Scholar CrossRef Search ADS   Park J., 1986. Synthetic seismograms from coupled free oscillations: Effects of lateral structure and rotation, J. geophys. Res.  91 6441– 6464. https://doi.org/10.1029/JB091iB06p06441 Google Scholar CrossRef Search ADS   Park J., 1987. Asymptotic coupled-mode expressions for multiplet amplitude anomalies and frequency shifts on an aspherical earth, Geophys. J. R. astr. Soc.  90 129– 169. https://doi.org/10.1111/j.1365-246X.1987.tb00679.x Google Scholar CrossRef Search ADS   Park J., 1990. The subspace projection method for constructing coupled-mode synthetic seismograms, Geophys. J. Int.  101 111– 123. https://doi.org/10.1111/j.1365-246X.1990.tb00761.x Google Scholar CrossRef Search ADS   Pekeris C., Jarosch H., 1958. The free oscillations of the earth, in Contributions in Geophysics in Honor of Beno Butenberg  eds Benioff H., Ewing M., Howell B., Press F., Pergamon Press. Resovsky J., Ritzwoller M., 1995. Constraining odd-degree Earth structure with coupled free oscillations, Geophys. Res. Lett.  22 2301– 2304. https://doi.org/10.1029/95GL01996 Google Scholar CrossRef Search ADS   Resovsky J., Ritzwoller M., 1998. New and refined constraints on three-dimensional Earth structure from normal modes below 3 mHz, J. geophys. Res.  103 783– 810. https://doi.org/10.1029/97JB02482 Google Scholar CrossRef Search ADS   Resovsky J., Ritzwoller M., 1999. Regularization uncertainty in density models estimated from normal mode data, Geophys. Res. Lett.  26 2319– 2322. https://doi.org/10.1029/1999GL900540 Google Scholar CrossRef Search ADS   Resovsky J., Trampert J., 2003. Using probabilistic seismic tomography to test mantle velocity-density relationships, Earth planet. Sci. Lett.  215, 121– 234. https://doi.org/10.1016/S0012-821X(03)00436-9 Google Scholar CrossRef Search ADS   Ritsema J., van Heijst H., Woodhouse J., 1999. Complex shear wave velocity structure imaged beneath Africa and Iceland, Science  286( 5446), 1925– 1931. https://doi.org/10.1126/science.286.5446.1925 Google Scholar CrossRef Search ADS PubMed  Ritsema J., van Heijst H., Woodhouse J., 2004. Global transition zone tomography, J. geophys. Res.  109 doi:10.1029/2003JB002610. https://doi.org/10.1029/2003JB002610 Ritsema J., Deuss A., van Heijst H., Woodhouse J., 2011. S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements, Geophys. J. Int.  184 1223– 1236. https://doi.org/10.1111/j.1365-246X.2010.04884.x Google Scholar CrossRef Search ADS   Rogister Y., Valette B., 2009. Influence of liquid core dynamics on rotational modes, Geophys. J. Int.  176 368– 388. https://doi.org/10.1111/j.1365-246X.2008.03996.x Google Scholar CrossRef Search ADS   Romanowicz B., 1987. Multiplet-multiplet coupling due to lateral heterogeneity: asymptotic effects on the amplitude and frequency of the Earth's normal modes, Geophys. J. R. astr. Soc.  90 75– 100. https://doi.org/10.1111/j.1365-246X.1987.tb00676.x Google Scholar CrossRef Search ADS   Tarantola A., Valette B., 1982. Generalized nonlinear inverse problems solved using the least squares criterion, Rev. Geophys. Space Phys.  20 219– 232. https://doi.org/10.1029/RG020i002p00219 Google Scholar CrossRef Search ADS   Trampert J., Deschamps F., Resovsky J., Yuen D., 2004. Probabilistic tomography maps chemical heterogeneities throughout the lower mantle, Science  306 853– 856. https://doi.org/10.1126/science.1101996 Google Scholar CrossRef Search ADS PubMed  Tromp J., Dahlen F., 1990. Free oscillations of a spherical anelastic earth, Geophys. J. Int.  103 707– 723. https://doi.org/10.1111/j.1365-246X.1990.tb05682.x Google Scholar CrossRef Search ADS   Um J., Dahlen F., 1992. Normal mode multiplet coupling on an aspherical, anelastic earth, Geophys. J. Int.  111 11– 34. https://doi.org/10.1111/j.1365-246X.1992.tb00551.x Google Scholar CrossRef Search ADS   Valentine A., Trampert J., 2012. Assessing the uncertainties on seismic source parameters: towards realistic error estimates for centroid-moment-tensor determinations, Phys. Earth. planet. Inter.  210–211 36– 49. https://doi.org/10.1016/j.pepi.2012.08.003 Google Scholar CrossRef Search ADS   Valentine A., Trampert J., 2016. The impact of approximations and arbitrary choices on geophysical images, Geophys. J. Int.  204 59– 73. Google Scholar CrossRef Search ADS   Woodhouse J., 1980. The coupling and attenuation of nearly resonant multiplets in the Earth's free oscillation spectrum, Geophys. J. R. astr. Soc.  61 261– 283. https://doi.org/10.1111/j.1365-246X.1980.tb04317.x Google Scholar CrossRef Search ADS   Woodhouse J., 1988. The calculation of the eigenfrequencies and eigenfunctions of the free oscillations of the earth and the sun, in Seismological Algorithms  pp. 321– 370, ed. Doornbos D., Academic Press. Woodhouse J., Dahlen F., 1978. The effect of a general aspherical perturbation on the free oscillations of the Earth, Geophys. J. R. astr. Soc.  53 335– 354. https://doi.org/10.1111/j.1365-246X.1978.tb03746.x Google Scholar CrossRef Search ADS   Woodhouse J., Deuss A., 2007. Earth's free oscillations, in Seismology and structure of the Earth  vol. 1 of Treatise on Geophysics, chap. 2, pp. 31– 65, eds Romanowicz B., Dziewoński A., Elsevier. Woodhouse J., Dziewonski A., 1984. Mapping the upper mantle: Three-dimensional modeling of Earth structure by inversion of seismic waveforms, J. geophys. Res.  89 5953– 5986. https://doi.org/10.1029/JB089iB07p05953 Google Scholar CrossRef Search ADS   Woodhouse J., Giardini D., 1985. Inversion for the splitting function of isolated low order normal mode multiplets, EOS, Trans. Am. geophys. Un.  66 300. Woodhouse J., Girnius T., 1982. Surface waves and free oscillations in a regionalized earth model, Geophys. J. R. astr. Soc.  68 653– 673. https://doi.org/10.1111/j.1365-246X.1982.tb04921.x Google Scholar CrossRef Search ADS   Woodhouse J., Giardini D., Li X.-D., 1986. Evidence for inner core anisotropy from free oscillations, Geophys. Res. Lett.  13 1549– 1552. https://doi.org/10.1029/GL013i013p01549 Google Scholar CrossRef Search ADS   Yang H.-Y., Tromp J., 2015. Synthetic free oscillation spectra: an appraisal of various mode-coupling methods, Geophys. J. Int.  203 1179– 1192. https://doi.org/10.1093/gji/ggv349 Google Scholar CrossRef Search ADS   APPENDIX: SYNTHETIC SEISMIC WAVEFORMS AND SPECTRA—THE NORMAL MODE FRAMEWORK A complete description of the theoretical basis for normal mode seismology lies well beyond the scope of this paper. Here, we aim to provide sufficient background to allow the reader to appreciate the motivations for—and distinctions between—self-, group- and full-coupling. In particular, we aim to highlight precisely where and why approximations need to be introduced within the various theories, providing insight into potential for systematic errors. Readers wishing for a more in-depth treatment are directed to a summary paper by Woodhouse & Deuss (2007), or the comprehensive account contained within Dahlen & Tromp (1998). Of course, our discussion here is very general: specific implementations of these techniques invariably adopt additional simplifications and restrictions, which may impact further upon accuracy. A1 The seismic wave equation The seismic wavefield within the Earth, s(x, t), can be described by the seismic wave equation. This has the general form   \begin{equation} \left(\mathcal {V}+\mathcal {W}\frac{\partial }{\partial t} +\rho \frac{\partial ^2}{\partial t^2}\right)\mathbf {s}=\mathbf {f} \end{equation} (A1)where $$\mathcal {V}$$ represents a linear integro-differential operator that depends upon the material properties of the Earth; $$\mathcal {W}$$ is an operator representing effects due to the Coriolis force; ρ = ρ(x) is the equilibrium density distribution within the Earth; and f = f(x, t) represents the force due to the seismic source. If viscoelastic effects are to be properly accounted for, the operator $$\mathcal {V}$$ is itself a function of time, and in this case $$\mathcal {V}\mathbf {s}$$ implies a convolution operation. Given specific realizations of the various quantities, we aim to solve eq. (A1) to obtain synthetic seismic observables. A general strategy for solving differential equations entails expanding the solution in terms of a known basis set, thereby reducing the equations to algebraic form. Although the solution remains the same irrespective of the basis chosen, the computational procedure required to obtain that solution can vary considerably. Thus, an appropriate choice of basis can be critical to developing efficient algorithms. In a seismological context, one strategy is to adopt a basis constructed on a finite-element–style mesh. This results in computational tools such as SPECFEM3D_GLOBE (e.g Komatitsch & Tromp 2002), where simulation costs scale with the length of time-series that must be calculated. Normal mode seismology takes a different approach: it relies upon basis functions that have global support, and the resulting algorithms have computational costs that depend largely on the desired frequency range of solution. Both approaches are formally complete and can provide faithful simulation of physical reality. However, we note that SPECFEM3D_GLOBE currently omits self-gravitation effects, as these are challenging to incorporate in a spectral element framework, and is therefore not suitable for modelling seismograms at the very low frequencies we consider in this paper (although efforts are being made to overcome this difficulty: see Gharti & Tromp 2017). A2 Spherical-earth eigenfunctions The particular family of basis functions we use are the eigenfunctions, or ‘normal modes’, of a spherically symmetric earth model. In this case, the material properties vary only with depth, and eq. (A1) simplifies considerably. The Coriolis term can be discarded, and the operator $$\mathcal {V}$$ takes a simpler form, which we denote by $$\bar{\mathcal {V}}$$. Similarly, we use $$\bar{\rho }$$ to highlight that the density distribution is also spherically symmetric. Neglecting the force term, and working in the frequency domain, it can then be shown that the differential equations admit a non-trivial solution, $$\boldsymbol {\mathcal {S}}_k(\mathbf {x})$$, at a discrete set of frequencies ωk. Thus, we have   \begin{equation} \left(\bar{\mathcal {V}}-\bar{\rho } \omega _k^2\right)\boldsymbol {\mathcal {S}}_k = 0 \end{equation} (A2)where $$\boldsymbol {\mathcal {S}}_k$$ represents the kth normal mode, and ωk is its corresponding eigenfrequency. This eigenvalue problem can be solved relatively straightforwardly (Pekeris & Jarosch 1958). There are an infinite number of $$(\omega _k, \boldsymbol {\mathcal {S}}_k)$$ pairs that satisfy eq. (A2), but only a finite set exist within any given frequency interval ω0 ≤ ωk ≤ ω1. Thus, it is feasible to compute a complete set of eigenfunctions within the given frequency range, and an algorithm described by Woodhouse (1988) provides a means for doing so. This method is implemented, for example, within the MINEOS software package described by Masters et al. (2011). It turns out (see e.g. Pekeris & Jarosch 1958) that distinct families of normal modes may be identified. Principally, two classes exist, known as ‘spheroidal’ and ‘toroidal’ modes; in addition, ‘radial’ modes are sometimes discussed as a separate category, although they are essentially a particular subset of the spheroidals. Furthermore, distinct sets of toroidal modes may be identified for each isolated solid region within the structural model, thus two groups of toroidal modes may be found for earth models such as PREM, where a solid inner core and mantle/crust region are separated by a fluid outer core. The modes take the spatial form of spherical harmonics laterally; their radial functions must be found by numerical integration. Thus, the index k on $${\boldsymbol {\mathcal {S}}}_k$$ may be recast as a set of four indices, k → (q, l, m, n). Here, q defines the ‘family’ of mode being considered; l and m are the angular degree and azimuthal order specifying the appropriate spherical harmonic Ylm(θ, ϕ); and n is the ‘overtone number’, which indexes the set of permissible radial functions. For present purposes, the main point of interest is that ωk depends only on q, l and n: in a spherically symmetric model, modes that differ only by their m-value share the same eigenfrequency. In the literature of normal mode seismology, it is common to see an individual (q, l, m, n) state described as a ‘singlet’; the (2l + 1) degenerate singlets corresponding to a given (q, l, n) triplet are then referred to as a ‘multiplet’. This triplet is often written in the form nql, with q being either S (in the case of spheroidal modes), or T (toroidals). For brevity, we will continue to make use of the index k within this discussion, unless there is reason to do otherwise. A2.1 Completeness and orthogonality The eigenfunctions of a spherically symmetric model can be shown to have some important properties, which stem from the self-adjointness of the operator $$\bar{ \mathcal {V}}$$. In particular, the infinite set of $${\boldsymbol {\mathcal {S}}}_k$$ form a complete basis for vector fields within the earth model, and they may be shown to be orthogonal, in the sense that   \begin{equation} \int _{V}\bar{\rho }(\mathbf {x}) \, {\boldsymbol {\mathcal {S}}}_j^*(\mathbf {x})\cdot {\boldsymbol {\mathcal {S}}}_k(\mathbf {x})\,\mathrm{d}^3\mathbf {x}=\delta _{jk} \end{equation} (A3)where δjk represents the Kronecker delta. Here, the integration volume V is the complete interior of the earth model, and an asterisk is used to denote complex conjugation. Strictly, these properties are only proven for entirely solid models, and we ignore any complications that might arise from the existence of ‘undertone modes’ in the Earth's fluid outer core. These are not thought to be problematic in a seismological context, as they exist at lower frequencies than we typically consider (e.g. Rogister & Valette 2009). A2.2 Seismograms in a spherically symmetric model Knowledge of the spherical-earth eigenfunctions allows straightforward synthesis of the wavefield expected within a spherically symmetric body. As shown by Gilbert (1970), and recounted in Woodhouse & Deuss (2007), the three-component wavefield can be reduced to the form   \begin{equation} {\bar{\mathbf {s}}}(\mathbf {x},t) = \sum _k \frac{1}{\omega _k^2}E_k (1-\exp \imath \omega _k t) \boldsymbol {\mathcal {S}}_k(\mathbf {x}) \end{equation} (A4)where the real part is understood and we only consider positive times. Ek represent the modal excitations, and may be obtained by calculating certain spatial integrals of the force term f(x, t). The process of computing seismograms by this formula is often referred to—for obvious reasons—as ‘mode summation’. In principle, the summation over k includes an infinite number of modes. However, we see that the dominant contribution of each mode to the seismogram lies at its eigenfrequency. Thus, if the seismogram is only required within a finite frequency band, only modes lying within that band need to be summed, and this makes the summation finite. In practice, we always work with band-limited seismic data, which is a necessary consequence of digital signal recording. If the maximum frequency present in our data set is ωmax, it is straightforward—and sufficient—to compute a mode catalogue containing all modes with ωk ≤ ωmax. A3 Seismograms in an aspherical earth model We now return to consider general, aspherical earth models, and the solution of eq. (A1). Invariably, the aspherical model will be specified as a 3-D structure superimposed upon a spherically symmetric reference model (often based on the model PREM, as described in Dziewonski & Anderson 1981). In what follows, we assume that the $$\boldsymbol {\mathcal {S}}_k$$ have been obtained using the relevant spherically symmetric reference structure. In addition, we implicitly assume that all interfaces within the aspherical model remain concentric spheres: in other words, no topography is present, either on the free surface or on internal discontinuities (such as the core–mantle boundary). To incorporate such topography, studies have traditionally relied on first-order boundary perturbation theory (Woodhouse & Dahlen 1978; Woodhouse 1980), which is inherently an approximation. Recently, an exact treatment of boundary perturbations has been developed by Al-Attar & Crawford (2016), although this theory has yet to be fully implemented. In any case, the effects of realistic-scale topography are not likely to impact the conclusions of this paper, and need not be considered further here. Given that the $$\boldsymbol {\mathcal {S}}_k$$ form a complete basis, we can express the wavefield s(x, t) in the form   \begin{equation} \mathbf {s}(\mathbf {x},t) = \sum _k u_k(t) {\boldsymbol {\mathcal {S}}}_k(\mathbf {x})\mathrm{,} \end{equation} (A5)where the real part is again understood. uk represent time-varying expansion coefficients. Strictly, this is only valid if the summation includes the entire and infinite set of eigenfunctions for the spherically symmetric model. However, we note the similarity with eq. (A4), and the fact that realistic aspherical models contain only relatively weak perturbations to the spherically symmetric reference model. It therefore seems reasonable to assert that modes with eigenfrequency far in excess of ωmax are unlikely to contribute within the frequency band of interest. We introduce a cut-off frequency, ωc ≥ ωmax, and work with the finite set of $$\boldsymbol {\mathcal {S}}_k$$ having ωk ≤ ωc. A3.1 Full-coupling In order to make use of this expansion, eq. (A5), we need to obtain expressions for the coefficients uk. Substituting into the wave equation, eq. (A1), we find   \begin{equation} \sum _k \left(u_k\mathcal {V}{\boldsymbol {\mathcal {S}}}_k+\dot{u}_k\mathcal {W}\boldsymbol {\mathcal {S}}_k+\rho \ddot{u}_k{\boldsymbol {\mathcal {S}}}_k \right)= \mathbf {f} \end{equation} (A6)and therefore, for any eigenfunction $$\boldsymbol {\mathcal {S}}_j$$, we obtain   \begin{eqnarray} &&{\sum _k\left(u_k \int _V \boldsymbol {\mathcal {S}}^*_j\cdot \mathcal {V}\boldsymbol {\mathcal {S}}_k\,\mathrm{d}^3\mathbf {x} + \dot{u}_k\int _{V}\boldsymbol {\mathcal {S}}^*_j\cdot \mathcal {W}\boldsymbol {\mathcal {S}}_k\,\mathrm{d}^3 \mathbf {x}\right.}\nonumber\\ &&{\quad+\,\,\left.\ddot{u}_k\int _{V}\rho \,\boldsymbol {\mathcal {S}}^*_j\cdot \boldsymbol {\mathcal {S}}_k\,\mathrm{d}^3\mathbf {x}\right) = \int _{V}\boldsymbol {\mathcal {S}}_j^* \cdot \mathbf {f}\,\mathrm{d}^3 \mathbf {x}.} \end{eqnarray} (A7)Since the sum only contains a finite number of terms corresponding to frequencies lower than ωc, this can equivalently be expressed as a matrix equation,   \begin{equation} \mathbf {Vu} + \mathbf {W\dot{u}}+\mathbf {P\ddot{u}}= \mathbf {q} \end{equation} (A8)where the elements of P, W, V and q are obtained by evaluating the appropriate integrals from eq. (A7). Taking the Fourier transform allows us to eliminate the derivatives, and introduce the matrix M(ω), such that   \begin{equation} \left[\mathbf {\tilde{V}}(\omega ) +\imath \omega \mathbf {W}-\omega ^2 \mathbf {P}\right]\mathbf {\tilde{u}} = \mathbf {M}(\omega ) \mathbf {\tilde{u}} = \mathbf {\tilde{q}} \end{equation} (A9)where the frequency-dependence of $$\mathbf {\tilde{V}}$$ arises from the inclusion of viscoelasticity within $$\mathcal {V}$$; strictly, eq. (A8) involves a convolution for the first term. The elements of M and $$\mathbf {\tilde{q}}$$ are readily evaluated, and thus obtaining $$\mathbf {\tilde{u}}$$ is a straightforward exercise in linear algebra. The ‘full-coupling’ approach involves solving eq. (A9) exactly. In reaching this point, we have introduced only one approximation: we have truncated the modal expansion at frequency ωc. To date, there has been little effort to investigate the errors introduced by this truncation, or assess how ωc should be chosen. Most of the existing studies that employ full-coupling assume it is sufficient to use ωc = ωmax, but there is little experimental or theoretical justification to support this. Thus, one goal of this paper is to explore how changes to ωc impact the seismograms produced, and to identify where the cut-off must be placed to ensure computations are essentially exact at frequencies below ωmax. We note that a number of ‘full-coupling’ studies introduce additional approximations within eq. (A9)—similar to those discussed in Appendix A4, below—designed to reduce computational complexity. It has been shown that these introduce only minor errors in results (Deuss & Woodhouse 2001). However, it is possible to avoid these simplifications, and the full-coupling calculations performed in this paper are based on solving eq. (A9) exactly. To achieve this, we adopt an iterative, preconditioned method described in Al-Attar et al. (2012), building on earlier work by Hara et al. (1993) and Deuss & Woodhouse (2004). A3.2 Self- and group-coupling Although solving eq. (A9) is conceptually straightforward, computational costs grow rapidly as ωc increases. Broadly speaking, the number of modes having ωk ≤ ωc (and thus, the dimension of the system in eq. (A9)) scales proportionally with $$\omega _c^3$$; the cost of obtaining $$\mathbf {\tilde{u}}$$ then grows as $$\omega _c^6$$. Computational costs rapidly become prohibitive, even with modern computational resources; for early studies, full-coupling was intractable. In consequence, it was necessary to identify strategies for simplifying the computational task. Looking again at eq. (A7), we can consider how the equations behave in the absence of any aspherical structure. In this case, all the integrals vanish unless j = k and thus, in the limit of a purely spherical model, the matrix M has non-zero elements only along its diagonal. This leads directly to our expression for seismograms in a spherically symmetric earth, eq. (A4), as is expected. However, it also motivates the assumption—borne out by empirical evidence—that for typical weakly aspherical models, M remains diagonally dominant. Elements lying far from the diagonal are, in general, relatively small—and in group- and self-coupling, they are therefore ignored. In eq. (A9), a general element of M is defined   \begin{equation} M_{jk}(\omega ) = \int _{V} \boldsymbol {\mathcal {S}}_j^*(\mathbf {x})\cdot \left[\mathcal {V} + \imath \omega \mathcal {W} - \omega ^2\rho (\mathbf {x})\right]\boldsymbol {\mathcal {S}}_k(\mathbf {x})\,\mathrm{d}^3 \mathbf {x}. \end{equation} (10a)In self-coupling, this is modified by introducing the assertion that   \begin{equation} M_{jk} = M_{(q_j,l_j,n_j)(q_k,l_k,n_k)}= 0 \,\,\textrm{unless}\,\, \left\lbrace \begin{array}{@{}l@{\quad }l@{}}q_j=q_k\\ l_j=l_k\\ n_j=n_k.\end{array}\right. \end{equation} (10b)In other words, self-coupling computes integrals using only the constituent singlets of each multiplet, and assumes that interactions between different multiplets can be entirely ignored. In group-coupling, the multiplets are assigned to groups. We introduce κ to denote a given multiplet, κ = {(q, l, m, n): m ∈ [−l, l]}. We then introduce groups, $$\mathcal {G}$$, containing one or more multiplets: $$\mathcal {G}_i = \lbrace \kappa _1^{(i)}, \kappa _2^{(i)},\ldots \rbrace$$. In general, all multiplets within a group will have similar eigenfrequency, although individual authors may differ over the particular groupings adopted. Group-coupling then replaces eq. (A10b) by the assertion that   \begin{equation} M_{jk} = 0 \,\,\textrm{unless}\,\, \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(q_j,l_j,m_j,n_j)\in \mathcal {G}\\ (q_k,l_k,m_k,n_k)\in \mathcal {G} \end{array}\right. \end{equation} (10c) for a single group, $$\mathcal {G}$$. Thus, some limited interactions are allowed between multiplets. However, (assuming a sensible ordering of the modes), both self- and group-coupling cause the matrix M to take on a block-diagonal form. Exploiting this structure allows a great reduction in the computational effort required to solve eq. (A9), at the expense, of course, of a deterioration in accuracy. A4 Normal modes in aspherical models Thus far, we have discussed the use of full-, group- and self-coupling to generate seismic spectra in aspherical earth models. These techniques make use of the eigenfunctions and eigenfrequencies of spherically symmetric models. However, it is also possible to compute the normal modes of general, aspherical models. One route to doing so is to adopt the same general strategy as we have already outlined: assume that the aspherical eigenfunction can be expressed in terms of the modes of a spherically symmetric reference model, and then solve for the necessary coefficients. This results in essentially the same analysis as in Appendix A3.1, except with time-invariant uk, and no force term. One may then regard eq. (A9) as a nonlinear eigenvalue problem, which can (in principle) be solved to yield the necessary expansion coefficients. The procedure is computationally challenging, limiting its range of application, especially in early studies. To simplify matters, we can make the same assumption as in self-coupling: that integrals of the form $$\int _V \boldsymbol {\mathcal {S}}_i^*(\mathbf {x})\cdot \mathcal {A}(\omega )\boldsymbol {\mathcal {S}}_j(\mathbf {x})\,\mathrm{d}^3 \mathbf {x}$$ are only non-zero when the singlets $$\boldsymbol {\mathcal {S}}_i$$ and $$\boldsymbol {\mathcal {S}}_j$$ correspond to the same multiplet. This allows the eigenvalue problem of eq. (A9) to be decomposed into a sequence of independent calculations, one for each multiplet. Furthermore, a perturbation-theory approach may be used to approximate the nonlinear eigenvalue problem as one in standard form. Given the orthogonality relationship, eq. (A3), we can write P = I + P1, where I is an identity matrix and P1 is the part arising only from the aspherical density structure (i.e. $$\rho -\bar{\rho }$$). We then assume that M(ω) can be expanded around a reference eigenfrequency of the multiplet, ω0, such that (cf. eq. A9)   \begin{equation} \mathbf {M}(\omega )\approx \mathbf {\tilde{V}}^{(\kappa \kappa )}(\omega _0)+\imath \omega _0\mathbf {W}^{(\kappa \kappa )} -\omega _0^2 \mathbf {P_1}^{(\kappa \kappa )} - \omega ^2\mathbf {I} \end{equation} (A11)where the superscript (κκ) is used to emphasize that each multiplet is being treated independently. This approximation relies on the assumption that |ω − ω0| is small, and that the earth model is dominated by spherical structure, allowing various additional terms involving P1, W and $$\partial \mathbf {\tilde{V}}/\partial \omega$$ can be neglected. It is common to introduce a ‘splitting matrix’, $$\mathbf {H} = \mathbf {\tilde{V}}(\omega _0)+\imath \omega _0\mathbf {W} -\omega _0^2 \mathbf {P_1}$$, allowing eq. (A9) to be written in the form   \begin{equation} \left(\mathbf {H}^{(\kappa \kappa )} - \omega ^2 \mathbf {I}\right)\mathbf {\tilde{u}}^{(\kappa \kappa )} = \mathbf {\tilde{q}}. \end{equation} (A12)Thus, the contribution of each multiplet to spectra—within the self-coupling approximation—is then found by diagonalizing the (2l + 1) × (2l + 1) matrix H. The (2l + 1) constituent singlets within the multiplet each contribute at a slightly different eigenfrequency, removing the degeneracy seen in spherically symmetric models. Thus, the spectral contribution of each multiplet becomes a set of (2l + 1) closely spaced—but distinct—peaks. This phenomenon is known as ‘splitting’. In practice, the resolution of spectral measurements may not allow individual singlets to be distinguished, but the overall shape of the peak corresponding to each multiplet depends upon the splitting, and thus upon the detailed properties of the Earth model. In the time domain, the contribution of the multiplet to seismograms can be written in the form   \begin{equation} \mathbf {s}^{(\kappa \kappa )}(\mathbf {x},t) = \mathbf {r}(\mathbf {x})\cdot \exp \left(\imath \mathbf {H}^{(\kappa \kappa )}t\right)\cdot \mathbf {v} \end{equation} (A13)where v and r are vectors that may be computed for the source and receiver, respectively; this result is due to Woodhouse & Girnius (1982). The discussion in this section has focused on self-coupling involving only individual isolated multiplets. Similar results can be obtained for group-coupling–assumptions: in effect, group-coupling can be regarded as self-coupling of a ‘super-multiplet’, allowing our analysis to translate straightforwardly. However, defining the reference frequency ω0 becomes more challenging, as the various multiplets within the group have different eigenfrequencies. Typically, some sort of average of these eigenfrequencies is used (for a recent discussion on specific averages, see Deuss & Woodhouse 2001). We note that the requirement for |ω − ω0| to remain small acts to restrict the size of the coupling group: only multiplets with similar eigenfrequency may be used. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial