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

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.


I N T RO D U C T I O N
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 1968Dahlen , 1969Woodhouse & Dahlen 1978;Woodhouse 1980;Woodhouse & Girnius 1982;Park 1986, Resolvability of Earth's 3-D Density 59 1987, 1990Romanowicz 1987;Tromp & Dahlen 1990;Lognonné 1991;Hara et al. 1991Hara et al. , 1993Um & Dahlen 1992;Deuss & Woodhouse 2001;Al-Attar et al. 2012;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 Earthlike 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  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. 1987Giardini et al. , 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.

T H E O R E T I C A L B A C KG RO U N D
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).

The forward problem
'Normal mode seismology' represents one framework for constructing a seismological forward model, particularly suited to lowfrequency, 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.
which describes the lateral form of 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 n q l (e.g. 0 S 2 ). The constituent eigenfunctions of each multiplet are then referred to as 'singlets'.

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) where S k is a vector field representing an eigenfunction of the reference model. There are an infinite number of eigenfunctions 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 Fouriertransformed expansion coefficientsũ k (ω) may be found by solving a system of linear equations (cf. eq. A9)

M(ω)ũ =q
( 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 vectorq is a representation of the seismic source (typically an earthquake).

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 .

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) where A(ω) is an integro-differential operator that depends upon frequency, and also upon the earth model. The elements tend to be relatively small unless S i and 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 This results in M taking a block-diagonal form, with each multiplet n q l 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 M i j (ω) = 0 unless S i and S j belong to the same group. (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;.

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.

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. (1987Giardini et al. ( , 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 κ = n q l 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.

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, κ = n q l and κ = n q l , 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 As discussed in Appendix A4, some additional assumptions allow us to introduce the splitting matrix, H, such that M (κκ ) ≈ H (κκ ) − ω 2 I. We then parametrize this by introducing a set of complex-valued 'splitting coefficients' c st such that where the γ are defined in terms of spherical harmonic triple products, 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' 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: η (κκ) , η (κ κ ) and η (κκ ) . 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 c st 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 with identical expressions for β and ρ. By deriving the exact expression for H (κκ ) mm (e.g. Woodhouse 1980), and then comparing this with eq. (5), it is possible to identify formulae for the kernels K such that 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).

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 1995Irving 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 dataprocessing 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 62 F. Akbarashrafi et al. 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.

T RU N C AT I O N E R RO R S I N F U L L -C O U P L I N G S E I S M O G R A M S
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 ?

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 (Swave 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'.

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 V s with respect to PREM, and where we added scaled compressional wave speed perturbations using δ ln V p = 0.5 δ ln V s and scaled density We refer to this as the 'reference' spectrum, s ref . 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. perturbations using δ ln ρ = 0.3 δ ln V s . 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.
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, s ref , and looked at amplitude spectra of the differences s (ωc=3 mHz) − s ref , etc., to provide an indication of the convergence of the fullcoupling 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 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 sρ − s 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 ) 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.
To quantify the truncation effects further, we defined a cumulative misfit for the entire synthetic data set 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 64 F. Akbarashrafi et al.  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 V p 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.  (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.

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 shortwavelength 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 V s 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 V p and ρ are scaled versions of V s , 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.
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.

Summary
The key conclusions we draw from these experiments may be summarized:  Figure 7. Self-coupling sensitivity kernels K (κκ) 00 (as in eq. 9) corresponding to the angular order s = 0 splitting coefficient for our selected modes. The solid red line denotes V p sensitivity, the solid black line V s 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.
(i) 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.
(ii) Accurate seismograms up to 3 mHz in earth-like models (having a red spectrum) require truncation at no less than 4 mHz.
(iii) Generally truncation errors are stronger at higher frequencies, as might be expected.
(iv) 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.

T H E I N T E R P R E TA B I L I T Y O F S P L I T T I N G F U N C T I O N M E A S U R E M E N T S
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. 2011Deuss et al. , 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. 2013Koelemeijer et al. , 2017 and density inferences (e.g. Ishii & Tromp 1999;Resovsky & Trampert 2003).
It is common for splitting function catalogues to publish uncertainty estimates on the recovered splitting coefficients and centrefrequencies. 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 selfor 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;, 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 (κκ) 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.

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 s 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 (κκ) 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: 0 S 2 , 0 S 5 , 0 S 6 , 0 S 7 , 1 S 4 , 1 S 7 , 1 S 8 , 1 S 9 , 1 S 10 , 2 S 6 , 2 S 12 , 2 S 13 , 5 S 3 . 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.

Observations from synthetic splitting functions
The coefficients of the splitting function measurements for mode 0 S 6 , a lower mantle mode, are depicted in Fig. 8. When selfcoupling 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.
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 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.
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 2 S 12 . 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.
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) where σ st are the uncertainties from Deuss et al. (2013) for the given mode and N represents the number of coefficients in the Reduced χ 2 misfit full−coupling seismograms self−coupling seismograms 2 S 6 2 S 12 2 S 13 5 S 3 0 S 2 0 S 5 0 S 6 0 S 7 1 S 4 1 S 7 1 S 8 1 S 9 1 S 10 Figure 10. 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).
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. 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 70 F. Akbarashrafi et al. Reduced χ 2 misfit δ in V s δ in V p δ in 2 S 6 2 S 12 2 S 13 5 S 3 0 S 2 0 S 5 0 S 6 0 S 7 1 S 4 1 S 7 1 S 8 1 S 9 1 S 10 Figure 11. Reduced χ 2 misfit defined by eq. (12) for model S20RTS, except with δ ln V s = 0 (black squares), δ ln V p = 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.
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 C tot = C th + C reg + C d , where C th is the theoretical uncertainty we wish to quantify, C reg is due to regularization in the regression and C d = σ 2 d 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 bootstrapping: randomly eliminating stations, and also whole events. This was done for a fixed regularization. Thus, these uncertainties therefore represent the data covariance, C 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, C reg = ασ 2 d 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 C th + C reg = βσ 2 d I. Hence β − α is a measure of C th expressed in units of the covariances in the Deuss catalogue. We see that C th and C 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 C tot = (1 + 1.3 + 2.0)σ 2 d I = 4.3σ 2 d 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 (eventstation 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 V s heterogeneity (having only lateral variations in V p and ρ), or no lateral V p structure (only V s 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 V s in the splitting functions is 2-3 orders of magnitude larger than the observational uncertainties for most modes, while that of V p 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 V s , around 4 for δ ln V p and around 1 for δ ln ρ.
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 V s 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.

Summary
Again, our results can be summarized as: (i) We confirm that there is a significant theoretical measurement error in currently published splitting functions.
(ii) There is no obvious bias, but this theoretical error can, to a first approximation, be assumed to be random noise.
(iii) There is also a significant contribution from regularization to the uncertainty in splitting functions.
(iv) Taking the total uncertainty into account, currently published splitting functions contain information on V s with a signal-to-noise ratio of about 100, information on V p with a signal-to-noise ratio of about 4, and information on ρ with a signal-to-noise ratio of about 1.
(v) Since the density signal is small in currently published splitting functions, advanced techniques should be employed to infer density structure.
(vi) The signal from V s 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 V s structure that may be derived these splitting function measurements.
These results are specific to the Deuss catalogue . 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.

C O N C L U S I O N S
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 selfor 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  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.

A C K N O W L E D G E M E N T S
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(FP/ -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.

A P P E N D I X : S Y N T H E T I C S E I S M I C WAV E F O R M S A N D S P E C T R A -T H E N O R M A L M O D E F R A M E W O R K
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 where V represents a linear integro-differential operator that depends upon the material properties of the Earth; 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 V is itself a function of time, and in this case Vs 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 timeseries 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 V takes a simpler form, which we denote byV. Similarly, we useρ 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, S k (x), at a discrete set of frequencies ω k . Thus, we have where 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 (ω k , 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 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 Y lm (θ , φ); 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 74 F. Akbarashrafi et al. 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 n q l , 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 selfadjointness of the operatorV. In particular, the infinite set of 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 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 where the real part is understood and we only consider positive times. E k 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 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 S k form a complete basis, we can express the wavefield s(x, t) in the form where the real part is again understood. u k 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 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 u k . Substituting into the wave equation, eq. (A1), we find k (u k VS k +u k WS k + ρü k S k ) = f (A6) and therefore, for any eigenfunction S j , we obtain 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, 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 Ṽ (ω) + ıωW − ω 2 P ũ = M(ω)ũ =q (A9) where the frequency-dependence ofṼ arises from the inclusion of viscoelasticity within V; strictly, eq. (A8) involves a convolution for the first term. The elements of M andq are readily evaluated, and thus obtainingũ 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