Can the Earth ’ s harmonic spectrum be derived directly from the stochastic inversion of global travel-time data ?

A set of seismic observations which all sample the same structure in the same way should have zero variance. This is naturally the case if all sources are in the same place, and the data are recorded by the same station. If sources and/or receivers are not in the same place, but close to one another, variance will generally be nonzero, but small. Variance might become large if the sampled region of the Earth contains heterogeneities whose spatial wavelength is comparable to the distances between sources and between receivers (and thus between the corresponding ray paths). The travel-time variance of a “bundle” of seismic rays thus reflects the degree of complexity of the sampled region of the medium. We apply this simple principle to real seismic databases, attempting to constrain the spherical harmonic spectrum of Earth’s structure without having to derive a tomographic model. This results in a reduction of the dimensionality of the solution space, and hence of computational costs. This approach allows to constrain the statistical properties, rather than exact geographic locations of structural features; knowing the statistics of Earth’s structure is most valuable for many fundamental geodynamic questions. We follow an earlier study by Gudmundsson et al. [1990] to find an approximate analytical relationship between averaged variance and harmonic spectrum; this allows us to determine the latter from a measurement of the former via a linear least squares inversion. Our analysis shows that the variance of ray bundles associated with large geographic extent of source/receiver bins is sensitive to low-degree spectral power, and vice-versa for small bins/high harmonic degrees. The method is accordingly ineffective at very low harmonic degrees, associated with an inherently limited number of source-receiver bins. We conduct a suite of inversions of both real and synthetic seismic data sets to evaluate the resolving power of our algorithm, and attempt to identify a range of harmonic degrees where the method is robust. Our results indicate that the resolution of the Earth’s spectrum afforded by the method presented here is inferior to that of classical tomography.


Introduction
After two decades of efforts to map the geographic distribution of mantle structure, the convergence between tomography and geodynamic models is only partial and limited to the larger scale lengths, while the small-scale components of Earth's structure are not well constrained [e.g., Becker andBoschi 2002, Bull et al. 2010].Whereas tomography remains the most widely employed tool to evaluate mantle structure, some authors have also implemented alternative methods, focusing on the statistical properties of mantle heterogeneity [e.g., Doornbos and Vlaar 1973, Haddon and Cleary 1974, Cormier 1999, Margerin and Nolet 2003, Garcia et al. 2009].
With this study we explore a "stochastic" approach alternative to tomography, introduced by Gudmundsson et al. [1990] (hereafter GDC90) and Davies et al. [1992] to constrain the overall strength of mantle heterogeneity as a function of depth, and estimate the variance of errors in teleseismic travel-time observations.The procedure of GDC90 allows to invert seismic observations to determine the depth-dependent sphericalharmonic spectrum of planetary structure, ignoring the geographic distribution of heterogeneity.This strategy is in principle useful because: (i) it involves a reduction of the dimensionality of the solution space: if, e.g., harmonic degrees up to 40 are considered, inverting for the harmonic spectrum rather than the 3-D structure of the Earth amounts to a two-order-ofmagnitude reduction of the number of dimensions in the solution space: this limits the non-uniqueness of the inverse problem, so that, particularly at high spherical harmonic degrees, the spectrum could in principle be constrained more robustly than it is now.(ii) The statistical properties, rather than exact geographic locations of structural features, are the piece of information that is most valuable for many fundamental geodynamic questions: general geodynamic models can reproduce only statistically the character of Earth structure.In this sense, the comparison of harmonic spectra obtained through modeling with those observed by seismology should provide valuable information on dynamic processes in the Earths interior [e.g.Bunge et al. 1996, Mégnin et al. 1997, Van Heck and Tackley 2008, Yoshida 2008, Foley and Becker 2009, Dziewonski et al. 2010, Nakagawa et al. 2010, Davies et al. 2012].
Importantly, our (and GDC90's) formulation requires that the stochastic process describing variations in seismic velocities in the Earth be Gaussian (i.e., velocity anomalies are normally distributed), isotropic (the correlation between two points depends only on the distance between them) over the entire mantle, and stationary at any given depth (the variance is the same for each point at that depth) (Section 3.1).
We apply our algorithm to two different global Earth-mapping problems: that of constraining global lateral variations in surface-wave phase velocity from teleseismic dispersion observations, and that of finding 3-D variations in P-wave velocity from a large traveltime database.Besides inverting real data, we evaluate the method's resolution with a suite of synthetic tests, aimed at identifying the range of harmonic degrees affected by the approximations required by the algorithm.

Stochastic formulation
A theory of wave propagation gives a mathematical relationship between relative anomalies in the properties of the Earth (e.g., the slowness p of a seismic phase) dp (r, i, {) (with r, i, { radius, colatitude and longitude, respectively) and anomalies dt in seismic traveltime.Neglecting non-linear effects, this relationship has the general form where V denotes the volume of the Earth, and the function K, dubbed sensitivity kernel (or partial derivative, Fréchet derivative), depends on the source-station geometry associated with the datum dt.Given phase and frequency, there exists one kernel per source-station couple.If a 1-D Earth is used as reference, the form of K depends only on epicentral distance.When ray theory is used to describe wave propagation, K is nonzero on the ray path (traced in the reference model), and zero everywhere else.If some form of finite-fre-quency theory is used, K becomes more complicated [e.g., Peter et al. 2007Peter et al. , 2009]].When the assumption of linearity is dropped, e.g. if we care about multiple-scattering, then no function K can be defined, and equation (1) ceases to be valid.Typically, equation ( 1) is used to set up an inverse problem with dp(r, i, {) as the unknown, and a set of observations of dt as data.
Equation ( 1) can be re-written for each observed value of dt, all of them with their corresponding kernel function K. dp is then expressed as a sum of unknown coefficients multiplied by some known "basis functions", e.g. with Y lm denoting the real scalar spherical harmonic of degree l and order m [e.g., Dahlen and Tromp 1998] and R n (r) some vertical basis function.The largest angular degree L and the total number of vertical functions N are selected depending on the resolution that one expects to achieve.Replacing (2) into (1) once per observation, we end up with a mixed-determined inverse problem with unknown coefficients A lmn , while all other quantities in (1) can be calculated.
The main idea of GDC90 and Davies et al. [1992] is to set up an inverse problem whose unknowns are not the coefficients A lmn , but the spectral power per unit area This is achieved by first defining an approximately equal-area grid spanning the Earth's surface.All dt observations associated with sources/receivers lying in the same pair of grid cells ("bins") are grouped in a "summary ray" [e.g., Morelli and Dziewonski 1987] or "ray bundle".To avoid possible biases related to grid geometry, this exercise is repeated four times, after as many rotations of the grid around the Earth's axis.The rotation angle coincides with the longitudinal extent of one of the equatorial equal-area grid cells, divided by five, so that after four rotations the entire longitudinal width of the grid cells is sampled.For each combination of horizontal grid size H, range of epicentral distance (distance "bin", identified by its mean value D) and range of source depth (depth bin, identified by its mean value Z), a value of the variance of dt is calculated, , , , , d , where k is the ray bundle index, from a total of n S ray bundles (taking into account also the bundles obtained after rotating the grid), n k is the number of actual rays collected in the k-th ray bundle, and mean k (dt) is the mean of all measurements of dt within that bundle.The factor n k in ( 4) is introduced so that v 2 is more strongly affected by summary rays formed by larger numbers of dt observations.For a given binning scheme (i.e., given values of H and Z), v 2 is a function of the epicentral distance D. Through (4), the numerical values of v 2 can be determined directly from the observations dt.For surface waves, we limit the summation over ray bundles only to those that include more than 10 rays, i.e. n k > 10.In addition, we consider only bins of (H, D, Z) with n S > 10.For body waves, we only take into account rays with travel-time |dt| < 4s and D < 100°, bundles with n k > 4 and bins with n S > 4, in analogy with GDC90.v 2 as defined by equation ( 4) can be thought of as the average, calculated over all ray bundles in the same (H, D, Z) bin, of the variance of dt calculated within each ray bundle.In practice, after introducing the operator (sum extended over all rays within a bundle) and the expected value operator (sum over all bundles in the 126 same (H, D, Z) bin), equation (4) takes the more compact form GDC90's theoretical treatment (pages 28 through 34) consists of showing that expression (5) can also be written as an integral function of the harmonic spectrum (3) of dp as a function of depth.That way, a linear inverse problem can be set up, whose unknowns are the coefficients Q ln themselves.In the following, we shall first rewrite all the theory for body waves in a more extensive way than GDC90 did (Section 3) and then reformulate it for surface waves (Section 4).

Formulation of the inverse problem for a 3-D Earth
Following GDC90, we take a stochastic approach, i.e. think of each ray bundle (for given (H, D, Z)) as a different realization of the same experiment.
After some algebra, equation ( 5) can be written The second term of the latter expression can be rewritten where A is the area of the grid cell of radius H and dt(x) is defined by equation (1).Inverting the order of the operators E and E C , and applying Fubini's theorem [e.g., Thomas and Finney 1996], we find Let us now consider the other term in equation ( 6), Since and again based on Fubini's theorem, If we define t = |x 2 − x 1 |, then dt(x 1 ) = lim t→0 dt(x 2 ), and Substituting (8) and ( 12) into (6), we find In Section 3.1 we shall show that, in the assumption of Gaussian, stationary and isotropic slowness perturbations dp, the expression E{dt(x 1 )dt(x 2 )}, which appears in both terms at the right-hand side of (13), can be written in a relatively simple form, function only of the distance t = |x 1 − x 2 |, and not of x 1 and x 2 themselves.In Section 3.2 we shall apply the raytheory approximation to further simplify the resulting expression.
In Sections 3.3 and 3.4 we use these results to (partly) solve analytically the double integral ∫ A ∫ A in equation ( 13), leading to a relatively simple expression for v 2 in terms of the harmonic spectrum of the Earth, Such equation constitutes the basis of GDC90's and our formulation of the inverse problem.

Relation between delay-time variance within ray bundles, and the harmonic spectrum of Earth's structure
We next show how E{dt(x 1 )dt(x 2 )} can be written in terms of the Earth's spectral coefficients Q l (r).Based on equation (1), with r i = (r i , i i , { i ) (i = 1,2) a position 3-vector defined within the Earth's volume V, and dr i the corresponding infinitesimal volume element.The kernel function K i is the one associated with the position 2-vector x i (i = 1,2) defined over the surface A, or the portion of Earth's surface swept by the ray bundle.
Equation ( 14) can be rewritten Let us focus on the term E{dp(r 1 )dp(r 2 )} at the right-hand side of this expression.Using equation (2), and denoting for simplicity We make at this point the important assumption, consistent with GDC90, that the fluctuations of dp are described, at any radius r within the mantle, by a Gaussian stochastic process, or in other words that, if we apply a shift Dr to the slowness perturbation map dp(r), the correlation between dp(r) and dp(r + Dr) quickly drops to zero with growing Dr.This property of dp(r) implies that E{dp(r 1 )dp(r 2 )} depends on the distance between r 1 and r 2 , and possibly their average radius r (the mantle's vertical coherence might vary with r), and a function f can be introduced such that with t the angular horizontal distance between r 1 and r 2 : We next write f (r,|r 1 − r 2 |, t) as a sum of Legendre polynomials P l (cos t) [e.g., Dahlen and Tromp 1998] with coefficients (2l + 1) c l (r, |r 2 − r 1 |)/4r, and The factor (2l + 1)/4r allows to simplify equation ( 19) after application of the addition theorem [e.g., Dahlen and Tromp 1998, eq. B.74], resulting in the expression Let us equate the alternative expressions for E{dp(r 1 )dp(r 2 )} found at the right-hand side of equations ( 16) and ( 21): It follows from ( 22) and the orthogonality of Y lm that and from (23) and, again, the orthogonality of Y lm that Following GDC90, we assume "some coherency in the harmonic pattern with depth", i.e. we assume that a function c exists such that A lm (r 1 )A lm (r 2 ) = A 2 lm (r) c(|r 2 − r 1 |)/(2l + 1), and equation ( 24) takes the form Note that we think of Q l as the unknown of an inverse problem, and make no distinction between Q l and its expected value E{Q l }.
We next substitute the expression (25) for c l into equation ( 19), and the resulting expression into (15), to find (26) a direct, linear relation between variance of dt within a ray bundle, and the Earth's harmonic spectrum

Simplification by application of ray theory, and the assumption that ray paths forming a ray bundle are parallel
As noted at the beginning of Section 2, in the raytheory approximation the velocity-kernel K(r) = 1 if r belongs to the ray path, and K(r) = 0 otherwise.Denoting ray 1 and ray 2 the ray paths corresponding respectively to the locations x 1 and x 2 within A, equation ( 26) can be rewritten Neglecting, at first, the effects of spherical geometry, GDC90 show in detail how the double integral ∫ ray 1 ∫ ray 2 can be reduced to a single integral along one reference ray.Their procedure requires the assumption that "all the rays have the same ray parameter and randomly distributed endpoints in the two grid cells defining the summary ray.This implies that the rays are approximately parallel and simply shifted horizontally with respect to each other" [GDC90, p. 32].Consider now a point r 1 on ray 1 .Let us call P its projection on ray 2 , and d the distance between r 1 and P (i.e., by the definition, the minimum distance between ray 1 and ray 2 ).Given a point r 2 on ray 2 , said s 2 the distance between P and r 2 along ray 2 , the distance x between r 1 and r 2 equals (see Figure 1).Then, for a function g(r 1 , r 2 ) that depends only on the distance x between r 1 and r 2 , The integral over x in equation ( 29) is particularly easy if g(x) = g 0 exp(−x 2 /a 2 ) for some a, implying Said x 1/2 the value of d such that g(x 1/2 ) = g 0 exp (x 2 1/2 /a 2 ) = 1/2 g 0 (i.e., x 1/2 is the "half-width" of the Gaussian g), it can be shown that x 1/2 = , and This result can be applied to the double integral at the right hand side of equation ( 27), assuming that its argument, a function of the distance t between points on ray 1 and ray 2 , be close to Gaussian.Then (32) Notice how the function c (|r1 − r2|) disappears between equation (26) and equation ( 32), once the assumption on the correlation between dp on the whole planet is assumed to be Gaussian.As already mentioned by GDC90, equation ( 32) holds approximately for many choices of an autocorrelation function, provided it does not have strong side lobes, i.e. the structure of the medium must not have a strong periodic component".Finally, in the assumption of parallel rays, the horizontal distance t has been systematically approximated with the generic distance d between ray paths, assumed constant along the ray paths themselves.

Writing the double surface integral as a single integral over distance
Recall the form of both the double surface integrals at the right-hand side of equation ( 13), i.e.
In Sections 3.1 and 3.2 we have rewritten the integrand E{dt(x 1 )dt(x 2 )}, expressing it in terms of the Earth's harmonic spectrum Q l (r) and reducing it to the simple form (32), function only of the constant distance between approximately parallel ray paths.Before using this result to set up an inverse problem, with Q l (r) as unknown and v 2 as datum, we reduce analytically the double surface integral∫ A ∫ A in (13) to a single, one-dimensional integral.

Cartesian case
ln cos x and d in equation ( 28).t is the horizontal distance between (i 1 , z 1 ) and (i 2 , z 2 ) introduced in equation ( 26). ( with A a circular surface of radius H, x 1 and x 2 Cartesian 2-vectors spanning A, and dx 1 , dx 2 the corresponding infinitesimal surface elements.Now, let the function f (x 1 , x 2 ) depend only on the distance t between the points x 1 and x 2 .We shall later make use of this property of f to simplify the double surface integral in equation (33).First, replace ∫ A dx 1 with an integral over the polar coordinates s (length) and | (angle), defined with respect to the centre of A. It follows that For each location (s, |), we must integrate again over all points x 2 within A. Let us replace the Cartesian coordinates x 2 with polar coordinates t (length) and ] (angle), defined with respect to the location x 1 , or (s, |).By definition, t then coincides with the distance between x 1 and x 2 .As we accordingly rewrite the integral ∫ A dx 2 in (34), we must specify the limits of integration in t and ]. t ranges between 0 and s + H.For each t, the interval of values of ] for which (t, ]) falls within A must be determined (and integrated over).If s the length z of the arc of ] to be integrated upon can be determined using the cosine rule, (see Figure 2 for a visual explanation of s, t and H).
Equation ( 34) now becomes where we have made use of the fact that f depends only on t, so the actual values taken by ] do not matter: only the length of the arc it spans does.For the same reason we can write f (s, |, t, ]) = f (t), and To further simplify the expression for I, it is convenient to change the order of integration over s and t.The first term at the right hand side of (37) We express the second term at the right hand side of (37) as the sum of two terms: one denoted I 2 , containing an integral over t between H − s and H, the other denoted I 3 , containing an integral over t between H and H + s.Changing the order of integration, we find (39) (40) If one now combines I = I 1 + I 2 + I 3 , equation (5) of GDC90 is reproduced.GDC90 further simplify the form of I, carrying out analytically the integrations over s.This is straightforward for the s-integral in I 1 , but more complicated for I 2 and I 3 .Those are solved via the formula with z = x 4 − b 2 x 2 + 2ax 2 + a 2 , valid for x > 0 and b > 0, leading to the final result of GDC90, 2 2 , cos cos ( (42)

Spherical case
To derive equation ( 42) we have treated the surface of the Earth, and of the area A spanned by a ray bundle, as flat.When their curvature is taken into account, equation ( 42) is replaced by equation ( 27) of GDC90 where and a, b 1 , b 2 are defined and

Relation between variance v 2 , and the harmonic spectrum of Earth's structure
We have shown in Section 3.1 that E{dt(x 1 )dt(x 2 )} is a function of the distance t = |x 2 − x 1 | between the two rays, assuming a parallel-ray approximation (i.e. the horizontal distance between two rays of the same ray bundle is approximately the same as the generic distance between the rays along their path).equations ( 33)-( 47) show that, since E[dt(x 1 )dt(x 2 )] = f (t), the double integral over surface is reduced to only one integral over distance by means of a weight function w(H, t), for both the cases of flat and spherical Earth.Replacing f (t) with E[dt(x 1 )dt(x 2 )] in equation (43), After substituting the term E[dt(x 1 )dt(x 2 )] with its explicit expression in (32), equation ( 48) becomes Let us now consider the first term at the right-hand side of equation ( 13).Using equation ( 32) with t = 0, we find where we have used the fact that P l (1) = 1, independent of l.Substituting ( 49) and ( 50) into ( 13), ( 51) where v 2 (H, D, Z) can be evaluated from the data, and at the right-hand side everything but Q l (r) is known.equation ( 51) constitutes the basis of the inverse problem solved by GDC90.

Projection to two dimensions
In a JWKB ray-theory description of surface-wave propagation [e.g., Ekström et al. 1997], where dv (i, {, ~) denotes lateral heterogeneities in surface-wave phase velocity v at angular frequency ~, and K is the corresponding sensitivity kernel.dv (i, {, ~) can naturally be rewritten as a linear combination of real spherical harmonics with constant A lm (no r-dependence).Here and in the following we shall neglect the dependence of dv, v, K and dt on ~; in practice, we shall always consider different frequencies separately.
We assume, as in the 3-D case, that the variations of dv (i, z) be Gaussian, stationary and isotropic.Equations (3)-( 13) remain then valid in the same form as above, provided that the radial basis-function index n, and the bin vertical extent Z, which are now meaning- (53) less, be removed.In a spherical reference frame, equation (15) can be rewritten In analogy with Section 4.2, we rewrite E{dv(i 1 ,{ 1 ) dv(i 2 ,{ 2 )} by either replacing dv 1,2 with their harmonic expansions, (which is the 2-D counterpart of equation ( 16)), or by invoking the statistical properties of dv (Gaussian, stationary, isotropic), which allow us to repeat steps ( 17) through ( 21) (after dropping the r-dependence of f ), resulting in with c l constant (no r-dependence).
After equating expressions ( 55) and ( 56) and repeating steps ( 23) through ( 25), we find the 2-D version of equation ( 26), with t denoting the angular distance between x 1 and x 2 .
In the ray-theory approximation, If, in analogy with Section 3.2, rays are treated as parallel, then a distance d between the two rays can be defined, and (see Figure 3 for a visual explanation of t, s 2 and d); equation (58) then becomes The t-integral in equation ( 60) is evaluated between d and r, because we shall only consider first-orbit surface-wave observations, hence t < r.Since Q l does not depend on t, equation ( 60) is rewritten (61) The t-dependent portion of equation ( 61) forms a convergent integral that can be solved analytically (Appendix A).We can swap summation and integration in (61), and define We show in Appendix A how to calculate c l (d).Substituting equation ( 62) into (61) and keeping in mind that, in the parallel-ray approximation, t = d, we are left with the compact expression

Relation between variance v 2 (H, D) and the harmonic spectrum of surface-wave phase velocity heterogeneity
As noted in Section 4.1, equation ( 13) is valid, in the same form, in both the body-wave and surfacewave cases.Making use of equations ( 61)-( 62), the second term at the right-hand side of (13) can be rewritten (64) According to equation ( 12), the first term at the right-hand side of (13) coincides with the limit of the second as t → 0. Let us focus first on the integrand at the right-hand side of ( 13): H dX 1 dX 2 .107) into (65), and the resulting expression, together with (64), into (13), we are left with The results of Section 3.3 apply, and the double integral over A can be rewritten as a single integral over distance, by means of a weight function w(H, t), so that equation ( 66) becomes (67) Since c l (t) does not depend on s, (68) which constitutes the basis of the surface-wave inverse problem that we shall solve in the following.

Solution of the inverse problem
We first address the 2-D inverse problem of Section 4, which is easier to implement, and later extend our formulation to the more complex 3-D problem of Section 3.

Discretization
We approximate the integral in equation ( 68 where D and F are the tensors defined by equations ( 74), (75), while Q is defined in Section 3.1.From equations ( 105) and ( 107), we notice that F n0 = 0 for all n.We then have no resolution at the harmonic degree l = 0, which will not be considered in all the following inversions.
It is common practice in the solution of linear problems to weight the input data vector D with a covariance data matrix C [e.g., Snedecor and Cochran 1980].
Assuming the values D n to be uncorrelated, the covariance matrix is diagonal, with diagonal entries coinciding with the weighted variance of the values of as defined in (4), that is [e.g., Bevington and Robinson 1992] Equation (76) then becomes (78) which can be solved in least-squares sense to find the spectrum nn and C − 1 -2 nn , respectively.

Least-squares solution and norm damping
We systematically discretize (H, D) so that N > L, where N is the largest possible value of n.Equation ( 78) is then an overdetermined problem which admits the least-squares solution [e.g., Trefethen and Bau 1997], where the superscript T denotes a transpose matrix, and we have made use of the fact that C − 1 -2 is diagonal and so Because seismic data are always polluted by measurement errors and their coverage is not uniform, the problem is ill-conditioned, i.e. the solution is not reliable unless equation ( 79) is regularized [e.g., Menke 1989].As a regularization constraint, we impose that the norm of the solution be minimum.The leastsquares formula (79) becomes ( 80) with m a regularization or "damping" parameter to be selected.
We solve equation ( 80) by means of Cholesky factorization of F T • C −1 • F [e.g., Press et al. 2001] (from now on LS) and of the non-negative least-squares algorithm (from now on NNLS) of Lawson and Hanson [1974].The solution Q, a spherical harmonic spectrum, is by definition positive, and NNLS guarantees that this constraint is satisfied.Cholesky factorization, on the other hand, has the advantage of being an exact method.
We apply the L-curve criterion [Hansen 1992] to select an adequate value of m: after defining the solu-tion norm (where the superscript m identifies the solution found with the corresponding value of the damping parameter), we divide it by the norm of the undamped solution o(0), thus defining the normalized norm (82) For each m, we also define the solution misfit and we can build the L-curve plotting the couples (o ~(m), g(m)).Our preferred value of m is the one corresponding to maximum curvature of the L-curve [Hansen and O'Leary 1993].

The 3-D problem: body-wave travel times and the depth-dependent spectrum of the mantle
We next discretize equation ( 51) to solve the original 3-D problem of GDC90.We first transform the ray integral over s to one over radius and find (84) with R bot the radius at the bottoming point of the ray and R ⊕ the Earth's radius.Then, if we choose a 1-D reference model and only consider a set of discrete values of (D, Z) (index i) and H (index j), it follows from equation (84) that where K i is the number of layers crossed by the ray associated to the i th bin and the length crossed by the ray through the k th layer.We -- establish a one-to-one correspondence between couples (i, j) and (k, l) and the values of two single indexes p and q respectively, so that in tensor notation equation ( 85) reads D = M • X, (86) with Since P 0 (x) = 1 for any value of x, it follows that M pq = 0 for l = 0, and, again, we cannot resolve and do not invert for the degree l = 0 spectral coefficient.
In analogy with Section 5.1.2,we solve equation ( 86) in a norm-damped least-squares sense, i.e.
(88) which we implement via NNLS, choosing m according to the L-curve criterion.As opposed to the treatment of Section 5.1.2above, and for consistency with GDC90, we do not weight the data through the covariance matrix in this case.

Application to global seismic databases
After calculating v 2 based on a set of surface-wave phase delays or body-wave travel-times, equations ( 80) and ( 88), respectively, provide the corresponding leastsquares solution for the surface-wave phase velocity or body-wave velocity spectrum.We implement (80) and (88) for two real global databases and compare the resulting harmonic spectra with those inferred, from the same data, based on tomography.Body-and surfacewave tomography maps are derived with the algorithm of Boschi and Dziewonski [1999].In this ray-theory/infinite-frequency approach, resolution is limited by the wavelength of inverted seismic and waves; the degree-40 spherical harmonic parameterization we utilize is within such limit, so that resolution is entirely determined by data "coverage" i.e. how well each H-D-Z bin is sampled by the data.

Surface-wave inversions
We apply our spectral inversion method to the fundamental-mode surface-wave dispersion database of Ekström et al. [1997], focusing for brevity on phase delays of Rayleigh waves at a period of 50 s (~65,000 observations) and of 100 s Love waves (~37,000 observations).
Figures 4a and 5a show the tomographic phase-ve-locity maps which we obtained from the same data, using a least-squares, ray-theory algorithm [e.g., Boschi and Dziewonski 1999], a spherical-harmonic parameterization up to degree 40, and the L-curve criterion to select regularization weight (roughness damping only).
The wavelength of the degree-40 zonal harmonic is ~1000 km, well above that of the longest-period waves considered here (100 s surface waves traveling at ~4 km/s, hence wavelength ~400 km) and thus the physical resolution limit of imaging.Figures 4a and 5a are in very good agreement with earlier results obtained from the same database, and so are the corresponding power spectra shown in Figures 4b and 5b; see in particular the maxima at degrees 2 and 5, corresponding to the ocean-continent signature, which, in this period range, e.g.Carannante and Boschi [2005] have found to be a robust feature, independent of the technique used to measure surface-wave phase dispersion.
In Figures 4b and 5b we compare the spectra found from tomography to those derived through our technique.The latter show a single maximum at degree 3, and generally more power than tomography at all harmonic degrees up to at least 10.The noise at high (>30) harmonic degrees is an effect of random noise in the data, as we have verified in synthetic tests on noisy data not shown here for brevity.The seemingly robust dominance of degrees 2 and 5 inferred from tomography is not found by the spectral inversion.The misfit is, nevertheless, low: 0.151 for the Rayleigh-, and 0.109 for the Love-wave inversion.The inversions of Figures 4b and  5b are regularized according to the L-curve criterion, but we have verified that changes in the weight of regularization do not particularly improve similarity to tomography results.

Compressional-velocity spectrum inversions
We implement equation ( 88) to determine the depth-dependent spectrum of mantle P-wave velocity from the data set of direct P-wave travel times of Antolik et al. [2001] and Antolik et al. [2003]; these are essentially International Seismological Centre P travel-time picks first selected by Engdahl et al. [1998] and then corrected by Antolik et al. [2003] for crustal heterogeneity (using the reference crustal model CRUST5.1 of Mooney et al. [1998]) after source relocation.As we are still focusing on the long-wavelength component of Earth's structure, data are collected in ~626,000 summary ray paths (2° bins) as described by Boschi and Dziewonski [1999].The "datum" v 2 , derived from P travel-times, is discretized as in the surface-wave case, and additionally binned according to source-depth; we select the following Z-bins: 0-30 km, 30-60 km, 60-100 km, 100-200 km, 200-450 km and 450-600 km, which coincide with those of GDC90.Values D p of v 2 are weighted according to the corresponding value of H p , so that the inversion is biased towards the low-degree spectrum, more difficult to constrain (repeating the inversion without this weighting results in a more unstable solution).Since Antolik et al. [2003] provides summary rays with 2° bins, no smaller value of H p are used in the calculation of v 2 .This does not pose a problem since a lateral resolution of ~2° roughly corresponds to harmonic degree l = 100, and we are concerned throughout this study with degrees l ≤ 40.
In analogy with Section 6.1, we invert the exact same database with the voxel-based mantle tomography algorithm of Boschi and Dziewonski [1999]; at each depth, we conduct a least-squares fit to find the degree-40 spherical-harmonic expansion that best approximates our voxel model.In Figure 6 we compare the tomography-based harmonic spectrum as a function of depth to that obtained from our spectral inversions.Differences between the two spectra are qualitatively similar to those between the surface-wave phase-velocity spectra of Sec-tion 6.1: at all mantle depths except for a few hundred kilometers below the transition zone, the tomography spectrum has a clear, well known [e.g., Becker andBoschi 2002, Dziewonski et al. 2010] maximum at degree 2; the "stochastic" spectrum much broader, especially at shallower depths, with the maximum centred at degree 3. The very broad spectrum at shallow depths was also observed by GDC90 (their Figure 17).Both tomographic and stochastic spectra show a change of character in the shallowest portion of the lower mantle, with the tomography maximum shifting from degree 2 to 1, and the stochastic one shifting from 2 to 5-6.Degrees 2 and 3 are again dominant at the bottom of the lower mantle.Again, the misfit is very low, g = 0.141.

Resolution analysis
Our method's failure to reproduce well established results of tomography is disappointing, but, as tomography is not error-free, it is not per se a proof of the method's ineffectiveness.We next evaluate directly the sensitivity and resolving power of the spectral method.74) describes the relationship between the unknown spectral coefficients Q l , and the value of v 2 associated with a (H, D) bin.On the basis of ( 74), let us introduce a sensitivity function (89) so that the actual kernel relating v 2 at H n with the spectral coefficient Q l coincides with the product D n K nl .

Surface-wave phase velocity spectrum resolution
We show in Figure 7 sensitivity K nl as a function of the size H n of geographic bins, and the harmonic degree l of the unknown spectral coefficients to be inverted for.Because in equation ( 74) the epicentral distance D n acts as a simple scaling factor, K nl alone fully describes the sensitivity of v 2 to Q l .As a general rule, we see from Figure 7 that different harmonic degrees are constrained by values of v 2 associated to systematically different geographic binning, i.e. large H n are needed to constrain the low-degree spectrum, while smaller H n serves to determine the higher-degree portion of the spectrum.It is immediately evident from Figure 7 that the averaged variance v 2 has little (but non-zero) sensitivity to harmonic degrees 1 and 2: it appears that the largest bin sizes employed here, 30° and 45°, are insufficient to provide sensitivity at degrees 1 and 2 comparable to the sensitivity of v 2 at degrees 3 and higher.At relatively high harmonic degrees, values of v 2 associated to a broad range of H n values are sen-sitive to the same spectral coefficient, so that fictitious coupling and trade-off between different harmonic degrees can be expected.

Resolution matrix
The resolution matrix associated with the inverse problem defined by equation ( 80 Figure 7. Sensitivity K nl defined by equation ( 89), as a function of harmonic degree l and bin size H.
our inversions: the closer R is to the identity matrix, the higher theresolution.Off-diagonal entries indicate "smearing" between the corresponding spectral coefficients.Diagonal values smaller than unit indicate that the spectral power associated with the corresponding harmonic degree might be underestimated [e.g., Boschi 2003].R depends on data coverage, but not on the data themselves.The value of the damping parameter m, on the other hand, has to be selected through the L-curve criterion as described in Section 5.1.2,and consequently depends on (the signal-to-noise ratio of ) the actual data that are inverted.
In Figure 8 we show R associated with the source/station list for the 50 s fundamental mode Rayleigh-wave data set of Ekström et al. verted in Section 6.1.We implement equation ( 90) via Cholesky factorization of F T • C −1 • F, using the value of m selected according to the L-curve criterion when applied to the inversion of real, 50 s Rayleigh-wave data (Section 6.1).R is far from the identity matrix, with values on the diagonal smaller than unit, suggesting a possible loss of amplitude.Figure 8 suggests that resolution at degree 3, where the diagonal entry of R is maximal, is higher than at degrees 1 and 2, confirming the poor sensitivity to low degrees seen in Section 7.1.1.Relatively large non-diagonal values, which we find in particular at high degree, indicate that neighbouring spectral coefficients will fictitiously map onto one another ("smearing"), again as anticipated in Section 7.1.1.

Spectral reconstruction of a random monochromatic model
In this and the next sections, we shall describe a suite of synthetic tests aimed at further quantifying the resolving power (or lack thereof ) of our algorithm.After defining a theoretical, "input" seismic velocity model, we calculate surface-wave travel-time delays dt in the ray theory approximation, i.e., we implement equation ( 52) integrating along the shortest great circle connecting source and receiver (no ray tracing).We then substitute the resulting synthetics into equation ( 4) to find the synthetic v 2 (H, D), and solve, again, the inverse problem (80).After the inversion, we compare the spectrum of the "input" model used to generate the synthetics with the one reconstructed by the inversion: their similarity is a measure of our algorithm's resolving power.
It should be noted that a synthetic v 2 (H, D) could alternatively be calculated by substituting the input spectrum Q l directly into (4): we verified that the two procedure yield consistent results (correlation r = 0.727 between the two resulting synthetic v 2 (H, D)).
Our first input model is constrained to have a simple monochromatic spectrum Q l = d l,9 : l,m coefficients are all 0 if l ≠ 9, while at l = 9 they are generated randomly.Based on the resulting model (Figure 9a), we calculate ~65,000 synthetic data (data set A), from the source/receiver distribution of 50 s Rayleigh-wave phase delays collected by Ekström et al. [1997] (Figure 9a).We do not add any noise to the data.We next compute the corresponding averaged v 2 , NNLS-invert it with our "spectral" algorithm and compare input and output spectra in Figure 9b.The maximum at degree 9 is roughly reconstructed, but smeared over a broad range of degrees, and with a drastic (order-of-magnitude) loss of amplitude.The performance of NNLS being so poor, we also look at LS-inversion (Section 5.1.2) results: lossof-amplitude is then not so severe, but strong aliasing occurs approximately between the dominant degree of the input model (l ~ 9) and its integer multiples.
We also apply the surface-wave tomography algorithm of Boschi and Dziewonski [1999], using a degree-40 spherical-harmonic parameterization as in Section 6.1, to invert data set A, and show in Figure 9b the resulting harmonic spectrum.We find that tomography also significantly underestimates spectral power, but perfectly reconstructs the monochromatic nature of the input model.The values of misfit g(m) associated with the spectral inversions in Figure 9b are g(m) = 0.458 (LS) and g(m) = 0.902 (NNLS); the misfit is large: compare, e.g., with the value of 0:072 found from the tomographic inversion of the same data with degree-40 harmonic parameterization.we generate 3,000,000 travel-time delays (data set B) associated with uniformly distributed sources and stations: this is a tremendous improvement in data coverage with respect to the 65,000, non-uniformly distributed observations of Section 7.1.3.In Figure 9b we compare the results of the subsequent inversions with the results of inverting the smaller data set A. We find the larger and more uniform coverage of data set B results in a better reconstruction of the maximum of the input spectrum, at l roughly between 8 and 10; in comparison with the results of tomography, however, resolution remains extremely poor.
We conclude that inadequacy of data coverage alone cannot explain the poor performance of the spectral inversion approach.7.1.5.Realistic phase-velocity spectrum We first derive a tomographic map of 150 s fundamental-mode Love-wave phase velocity (Figure 10a), this time inverting ~16,000 dispersion observations from the database of Trampert andWoodhouse [1996].The tomography algorithm is the same as in Section 6.1, with sphericalharmonic parameterization up to degree 40, regularized through simple norm-damping and according to the L-curve criterion.
From the phase-velocity map of Figure 10a we calculate a set of synthetic phase delays, associated with the same source/receiver distribution of the original, real data set.We then compute the corresponding v 2 (H,D) and invert it with our algorithm.We compare in Figure 10b the resulting output spectrum with that of the input (tomography) model.The two peaks at l = 2 and l = 5 are merged into a single maximum of the output spectrum at l = 3: this reminds one of the discrepancy between stochastic and tomography spectra obtained from real data in Sections 6.1 and 6.2.As in Section 7.1.3the misfit is high (g(m) = 0.540 against ~0.05 achieved by tomography).In summary, application of our algorithm to a realistic, though noise-free, data set confirms the negative results of Sections 7.1.3and 7.1.3.

Body-wave velocity spectrum resolution
7.2.1.Sensitivity of v 2 to Q l Equation (29) of GDC90 describes the sensitivity K 3D nl of v 2 , for a given grid-size H, to the spectral coefficient Q l of mantle heterogeneity, i.e. (91) The kernel K 3D nl is the 3-D counterpart of K nl as defined by equation ( 89) above.Discretizing as explained in Section 5. 1.1, (92) We plot K 3D nl in Figure 11 as a function of bin size H n and harmonic degree l. Figure 11 indicates that the P-wave database is, like the surface-wave one, only marginally sensitive to degree-1 and -2 structure, independent of depth in the mantle.Compared to the surface-wave case (Figure 7), values of v 2 associated with large H n are sensitive to structure at relatively large harmonic degrees (e.g. for H n = 30°, K 3D nl at l = 5 is about as large as K 3D nl at l = 40.Because sensitivity of v 2 is high over a large range of harmonic degree, alias/smearing can be expected in spectral inversions.Note that Figure 11, including the high sensitivity of v 2 at large H to high degree Q l , is in good qualitative agreement with Figure 10b of GDC90.

Resolution matrix
The resolution matrix associated with equation ( 88) is (93) Again, R does not depend directly on the data (though it does depend on their geographical coverage).It depends on the data indirectly through m, which is selected according to the L-curve criterion (Section 5.1.2),and is affected by the signal-to-noise ratio of the actual observations (larger noise requires stronger damping).We parameterize the statistics of mantle structure in terms of 10 harmonic spectra, 1 ≤ l ≤ 40, each associated to one of 10 equal-thickness layers.The matrix is numbered so that the first 40 indexes are associated to the 40 values of Q l (l = 1, ..., 40) at the top layer the following 40 to the second shallowest layer and so on down to the bottom of the mantle.This results in R being approximately block-diagonal, with as many blocks as there are layers in our vertical parameterization.Each block on the diagonal corresponds to resolution of an individual layer, and smearing within the layer.Off-diagonal blocks correspond to smearing between the same or different harmonic degrees in different layers.
We show in Figure 12a R calculated from the ~626,000 source-station couples of Antolik et al. [2003] (Section 6.2).Entries within the diagonal blocks are systematically much larger than throughout the rest of the matrix.This indicates that the coupling between spectral coefficients at different degrees and depths is limited, and "vertical" resolution acceptable.If lateral resolution were high, i.e. the coupling between different harmonic were low, diagonal blocks would be closer to diagonal, in analogy with Figure 8.This is not the case: many off-diagonal entries within each diagonal block are comparable to diagonal ones, indicating that Q l is poorly resolved: smearing/fictitious coupling between harmonic coefficients within a layer occurs (Figure 12b).Resolution is poor independent of l, i.e. unlike tomography, the spectral method does not a better job of resolving low harmonic degrees than it does with high degrees.As anticipated, it is inherently difficult to project information on very-low-degree mantle structure into the function v 2 .

Spectral reconstruction of a random monochromatic model
We employ a vertically homogeneous input model, with a pattern of P-velocity variations coincident, at all .

R M M
/ DELLA MORA AND BOSCHI depths, with the model of Section 7.1.3.We generate ~626,000 synthetic P travel-time delays making use, again, of the source-station couples of Antolik et al. [2003] and adopting a linear relation between model anomalies and data [e.g., Boschi and Dziewonski 1999].We calculate v 2 (H, D, Z) via equation ( 4), with M now defined by equation ( 87).Consistently with Section 7.1.3,we verify that this is approximately equivalent to substituting Q in equation ( 86) with the input-model spectrum.
The results of least-squares solving the inverse problem (88) are shown in Figure 13, and are characterized by the same resolution problems encountered in the 2-D case.The impulsive nature of the spectrum is not reproduced at any depth.The smearing around the main peak is comparable with that of Figure 9, confirming that the spectral method cannot effectively discriminate between individual harmonic degrees.The achieved misfit g(m) = 0.267 is good, which, together with the strong discrepancy between input and output spectra, indicates that the sensitivity of v 2 (H, D, Z) to the Earth's spectrum is severely limited.

Discussion and conclusions
With this study we attempted to devise algorithm to estimate the complexity of planetary structure, defined in a spherical harmonic parameterization, directly from a global set of seismic observations.Our procedure relies on the assumption that Earth heterogeneity be adequately described as a Gaussian, stationary and isotropic stochastic process.It is based on the idea that, if one subdivides the Earth's volume into a set of regions of a given size (each interpreted as an independent realization of the same experiment), the mean variance of dp within regions is related to the spectral power of spherical harmonics of wavelength comparable to the size of the region.The number of possible independent subvolumes into which the Earth can be subdivided decreases with the harmonic degree,  so that the low-degree spectrum is inherently undersampled and presumably poorly resolved.This problem has an analogous in the estimation of cosmic properties at scales close to that of the observable universe: in cosmology, there is a large uncertainty on these quantities, based on the idea that it is possible to have many independent observations (and therefore perform a statistical analysis on them) only for small-scale properties of the observable universe, whereas this does not hold for structures that are comparable with its size [e.g., Somerville et al. 2004].There is nevertheless no a-priori reason to exclude that the intermediate-and higher-degree spectrum of the mantle can be constrained by very large, recent seismic databases.
While no measure and/or inversion method can provide a ground-truth observation of the Earth's spectrum, some of its properties are by now clearly well constrained and generally accepted.In the uppermost mantle, the robustness of the Earth's spectrum up to degree ~12 is argued for by, e.g., Carannante and Boschi [2005], who found highly correlated results from completely independent databases; for relatively shortperiod (~30 s) surface waves, for example, spectral peaks at degrees 2 and 5 are found independent of the observation and inversion techniques.The lower-mantle spectrum has been analysed by Becker and Boschi [2002], who find common spectral features from a wide variety of P-and S-velocity tomography models.Particularly robust is the dominance of degree 2 at most mantle depths.The conclusions of Becker and Boschi [2002] on at least the long-wavelength portion of the Earth's spectrum have been confirmed by more recent global-tomography studies.None of those features is reproduced by the surface-and body-wave spectral inversions illustrated here, which systematically result in smoother spectra without sharp maxima at any, low or high harmonic degree.While the results of tomography cannot be taken as ground truth, even at the longest wavelengths, the disagreement with such well established features suggests that the spectral approach might simply lack the resolution needed to properly extract the Earth's spectrum from seismic data.
The ineffectiveness of the spectral method is confirmed by the resolution analysis of Section 7. We have designed ad hoc synthetic experiments to try and determine which particular simplification in GDC90's and our formulation could limit resolution so severely.One important assumption is that source/receiver coverage be sufficiently uniform to sample Earth's structure at all scales.This is probably not the case in the real world, since sources are essentially limited to plate boundaries, and receivers to continents and ocean islands.In Section 7.1.4we illustrate the results of a synthetic test in-volving millions of data from an uniform source-station distribution.Even in such an idealized scenario, the a-priori spectrum is far from being recovered (Figure 9b).
We explored several other possible reasons for the failure, or lack of resolution of our inversions.The regularized linear inversion procedure per se can generate artifacts, described by the resolution matrices shown in Figures 8 and 12: this can in principle be avoided through a nonlinear inversion procedure, but we have verified that nonlinear inversion (genetic algorithm) practically does not improve the resolution of our method, as documented in detail by Della Mora [2012]; the nonlinear approach also allowed Della Mora [2012] to drop some of the approximations required here; yet, resolution remained equally poor [Della Mora 2012, § 2.8].The very low sensitivity of v 2 to the low-degree spectrum, inherent to our problem as mentioned above, is confirmed quantitavely by Figures 7 and 11.At intermediate degree, higher resolution could in principle be possible, but aliasing of unresolved, lower-degree signal is presumably an issue, resulting in poor resolution at all harmonic degrees.We speculate that aliasing could be limited, and intermediate-wavelength complexity more robustly constrained, with a choice of basis functions different from spherical harmonics, but this question is better left to future research.$.   112) and the numerical one of equation ( 108). (105)

A2. Numerical implementation -Validation of our analytical expression for c l
Before the inverse problem ( 80) is solved, we must calculate the numerical values of the matrix entries F nl .This involves the calculation of the function c l according to its analytically derived expression (105).
We validate our analytical integration of equation (62), carried out in Section A, by comparing its result ( 103) with the values of c l found from equation (94) (l = 0) and by direct numerical integration of equation ( 95).We do not integrate equation ( 62) numerically because it is singular at t = d.Substituting equation ( 96) in ( 95), (108) In our approach, we substitute in equation ( 105 The advantage of our approach in equation ( 112) with respect to the integration of equation ( 108) is that the term oscillates much more intensively then sin [(l − 2k) t] and, because of this, it gives less numerical problems to integrate the latter function rather than the former.
The result of this comparison is summarized in Figure 14: it can be seen that the error has the maximum values in the same points of the analytical expression of equations ( 94)-( 95), and looking at the colour scales of the two plots it is evident that the error is approximately six orders of magnitude smaller than the exact value.

Figure 2 .
Figure 2. Sketch of variables s, t and H in equation (35).A is the surface where the integrals in equation (33) are done, C is the centre of A and H and s are the same as in equation (35).
) with a discrete summation, to find with We next discretize values of cell size H and angular distance D, and equation (69) takes the form (71) While the discretization of D is straightforward (we simply subdivide the domain D = 0° < D < 180°i nto 2° intervals), that of H is more problematic: we want each of the four sectors of the Earth's surface defined by the equator and the Greenwich meridian to contain an integer number of grid squares.This results in the constraint(72)    where mod[a,b]  is the remainder of the division of a by b, with a ∈ N and b ∈ N, and int(x), with x ∈ N, is the largest integer not greater than x.In practice, we employ 39 indexed values of H in the range 3° ≤ H ≤ 45°, ordered by increasing value.Values are closely spaced (≤ 1°) up to H = 15°, after which only H = 22.5°, 30°, 45° is possible.Values of D are regularly discretized from D = 1° to D = 179° with a sampling of 2°.A one-to-one correspondence is established between couples (i, j) and the values of a single index n,

Figure 4 .
Figure 4. (a) Tomography map of 50 s fundamental-mode Rayleigh-wave phase-velocity, with superimposed sources (red squares) and stations (green circles) of the inverted database.(b) Result of spectral inversion (blue curve), compared with the spectrum inferred from the tomography map (black).

Figure 5 .
Figure 5. Same as Figure 4, but for the 100 s fundamental-mode Love-wave data set.

Figure 6 .
Figure 6.Spherical-harmonic spectrum as a function of depth (a) from global tomography applied to a large P-wave database, and (b) from the stochastic inversion of the exact same data.Both spectra are independently normalized by the maximum power at each depth.

Figure 8 .
Figure 8. Resolution matrix from equation (90) for the LS solution for a dataset of Rayleigh waves at 50 s from Ekström et al. [1997].
Figure 9. (a) Monochromatic random input model with Q l = d l,9 .Red squares represent sources and green circles represent stations.(b) Spectra of the input (black curve) model of panel 'a', and of the output models resulting from LS (red) and NNLS (blue) inversions of synthetic data sets A (solid) and B (dotted) generated from the model of panel 'a'.Note the different scales for input and output models.

Figure 10 .
Figure 10.(a) Tomography model of 150 s fundamental-mode Love-wave phase velocity, with superimposed source (red squares) and station (green circles) locations.(b) input (from the model of panel 'a') and output harmonic spectra, associated with the synthetic test described in Section 7.1.5.

Figure 11 .
Figure 11.Same as Figure 7, but with the 3-D sensitivity function of equation (92).
Figure 12.(a) 3-D, P-wave resolution matrix defined in Section 7.2.2, computed based on the database of Antolik et al. [2003].(b) The block framed in black in panel 'a', corresponding to the third layer from the top.

Figure 13 .
Figure 13.Depth-dependent spectrum resulting from the spectral inversion of our synthetic P arrival-time data set.The data set includes ~626,000 synthetic observations based on a vertically homogeneous input model, with the pattern of P-velocity variations coincident with the random degree-9 model of Section 7.1.3.In case of perfect resolution, Q l should be zero at all harmonic degrees except for l = 9, at all depths.
(102)  with respect to t and substituting the result into (95) (103)In the following, we shall compact the notation by defining (104)Combining equations (94) and (103) we find the following expression for the integral (62):Finally, equation (66) requires the calculation of lim t→0 c l (t).From the definition of c l (t), (106) Using the expression (102) for P l (cos t), this re-

Figure 14 .
Figure 14.Validation of our analytical expression (112) for c l ; (left panel), values of c l from equation (112); right panel, the difference between c l from the analytical evaluation of equation (112) and the numerical one of equation (108).
) the definition of h l−2k (d) of equation 104, obtaining Finally we calculate numerically h l−2k (d) with the approximated formula where so that the resulting expression of c l (d) is