Jérôme Mariette
email: jerome.mariette@toulouse.inra.fr

Nathalie Villa-Vialaneix
email: nathalie.villa@toulouse.inra.fr

Aggregating Self-Organizing Maps with Topology Preservation

Keywords: Self-Organizing Maps, aggregation, topology preservation

In the online version of Self-Organizing Maps, the results obtained from different instances of the algorithm can be rather different. In this paper, we explore a novel approach which aggregates several results of the SOM algorithm to increase their quality and reduce the variability of the results. This approach uses the variability of the algorithm that is due to different initialization states. We use simulations to show that our result is efficient to improve the performance of a single SOM algorithm and to decrease the variability of the final solution. Comparison with existing methods for bagging SOMs also show competitive results.

Introduction

Self-Organizing Maps (SOM), [START_REF] Kohonen | Self-Organizing Maps, 3rd Edition[END_REF] have been shown to be powerful methods for analyzing high dimensional and complex data (see, for instance, [START_REF] Kohonen | MATLAB Implementations and Applications of the Self-Organizing Map[END_REF] for applications of the method to many different areas). However, the method suffers from its lack of good convergence properties. In its original version, the theoretical convergence of the algorithm has only be proved in very limited cases [START_REF] Cottrell | Theoretical aspects of the SOM algorithm[END_REF] and even in the modified version in which the training of the SOM is expressed as an energy minimization problem [START_REF] Heskes | Energy functions for self-organizing maps[END_REF], different runs of the algorithm give different results, that can be very dependent on the initialization. This problem is even more critical when the data set to be analyzed is complex or high dimensional.

This paper addresses the issue of aggregating several results of the SOM algorithm, all obtained on the same data set. Several attempts to combine SOMs while preserving their topological properties have been proposed in the literature [START_REF] Petrakieva | Bagging and bumping self organising maps[END_REF][START_REF] Saavedra | Fusion of self organizing maps[END_REF][START_REF] Vrusias | Distributing SOM ensemble training using grid middleware[END_REF][START_REF] Baruque | Fusion methods for unsupervised learning ensembles[END_REF][START_REF] Pasa | Fusion of Kohonen maps ranked by cluster validity indexes[END_REF]. In this paper, we present a novel method to combine several SOMs while preserving their topology. The proposed method combines several ideas taken from the different methods and allows to explore initialization states. It is both simple and efficient. We present a full comparison of the different options to aggregate the results of different SOMs and discuss the most relevant choices. Finally, we show that our approach is a competitive alternative to the existing methods on real data applications.

The remainder of the paper is organized as follows: in Section 2, an overview of aggregation methods for SOMs is presented. In Section 3, the proposed method is described. Finally, Section 4 presents experimental results and comparisons.

An overview of aggregation methods for SOMs

Suppose that B results of the SOM algorithm are given for the items (x i ) i=1,...,n , (M b ) b=1,...,B . Each of these results, M b is well defined by its set of prototypes (p b u ) u=1,...,U and comes with an associated clustering function

φ b : x ∈ R d → arg min u=1,...,U x -p b u 2 .
For the b-th SOM, the clusters will be denoted by (C b u ) u=1,...,U , where

C b u = x i : φ b (x i ) = u .
The purpose is to build a fused or a merged map, M * , with prototypes (p * u ) u=1,...,U and a clustering function φ * which improves and summarizes the B maps into a unique consensual map. Note that all SOMs have been trained from the same data (x i ) i=1,...,n or from a subset (e.g., a bootstrap sample) of this data set. They can also have been trained from different descriptors of the observations (e.g., from different sets of variables observed on the same items): in this case, the fused map thus corresponds to a map integrating the different descriptors. However, for the sake of simplicity, we will restrict our description and simulation to the first case (same observations, or eventually, bootstrap samples from the same observations and same descriptors).

As already explained in [START_REF] Petrakieva | Bagging and bumping self organising maps[END_REF] in the context of a one-dimensional grid, there is no ground truth for cluster labelling in the unsupervised framework. A first strategy to overcome this issue is to perform a re-labelling of the clusters based on the clustering only: [START_REF] Saavedra | Fusion of self organizing maps[END_REF] merge together the clusters of different maps with a majority vote scheme. A "fused" prototype is defined as the centroid of the grouped cluster prototypes over b = 1, . . . , B and a topology is deduced posterior to the definition of the clusters. Another approach that uses the different maps in an indirect way is described in [START_REF] Mariette | Bagged kernel SOM[END_REF]: in this paper, we proposed to use a subset of (x i ) i , using the most representative observations of the set of B maps, to train a final SOM from a simpler and more robust data set. This method is well suited to handle very large data sets. However, both approaches do not necessary produce a map with a topology similar to the B merged SOMs and make use of only a small part of the information provided by the B learned SOMs.

Several attempts to explicitly take advantage of the prior (common) structure of the maps have been proposed in the literature. A first method consists in constraining the B SOMs to be as similar as possible by a common initialization. This initialization can be derived, for instance, from a PCA of (x i ) i . Then, the different maps are fused by averaging the prototypes of the clusters situated at the same position the B SOMs [START_REF] Vrusias | Distributing SOM ensemble training using grid middleware[END_REF] or by using a majority vote scheme to classify the observations [START_REF] Petrakieva | Bagging and bumping self organising maps[END_REF]. Alternatively, [START_REF] Petrakieva | Bagging and bumping self organising maps[END_REF][START_REF] Baruque | Fusion methods for unsupervised learning ensembles[END_REF] also propose to make the B SOMs similar by initializing the b-th SOM with the final prototypes of the previous one. [START_REF] Baruque | Fusion methods for unsupervised learning ensembles[END_REF] improves this approach by weighting the averaging of the prototypes by a cluster quality index. Similarly, [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF] uses a similar strategy to handle streaming or large data sets, splitting the data into several patches that are sequentially processed by a different SOM algorithm initialized with the result of the previous one. However, these methods do not allow to explore the possibilities of different initializations, which can be an issue in SOM. Moreover, a sequential initialization of the B SOMs prevents from training them in parallel, which can be an important issue if B is large: using a large B is advised for stabilizing the result of the algorithm.

Another approach to preserve the topology property of the map is to align the different maps on one of them, which serves as a reference for the topology: in [START_REF] Georgakis | An ensemble of som networks for document organization and retrieval akrr[END_REF], the map is chosen arbitrarily, and the other maps are fused sequentially to this first one, averaging the prototypes (p b u ) u of the current map to the closest prototypes of the current fused map (p * u ) u . To leverage the problem of the choice of the map that is used to align the other maps, [START_REF] Pasa | Fusion of Kohonen maps ranked by cluster validity indexes[END_REF] proposes to choose a reference map that is the best one according to a given clustering quality criterion. However, this method makes the result strongly dependent on the choice of the first map because only its topology is used, whereas the topologies of the next maps are not utilized as such.

Description of the optimal transformation method

It is well known that the quality of the SOM strongly depends on its initialization. Given different maps obtained from different (random) initializations, we propose to find the "best" transformation that can be used to obtain two comparable results between two distinct maps. The optimal one-to-one transformation between prototypes in general might be difficult to define so we restrict ourselves to transformations that strictly preserve the topology of the map, i.e. the set of linear isometric transformations (rotation and/or symmetry). To do so, only square maps with m rows and columns are considered (i.e., using the notations introduced in the previous section, U = m 2 ): in these maps, the clusters are supposed to be positioned on a 2D grid at coordinates

{(k 1 , k 2 )} k1, k2=1,...,m .
Then, T denotes the set of all transformations, T : R 2 → R 2 , that let the map globally invariant: more precisely, T is composed of the set of rotations {r θ } θ∈{0,π/2,π,3π/2} and of the transformations {r θ • s} θ , with s the symmetry with respect to the axis passing by the points m+1 2 , 0 and m+1 2 , m . For a given map M with prototypes (p u ) u and a given T ∈ T , the transformed map T (M) is the map in which the unit u, with coordinates (k u 1 , k u 2 ) in N 2 , has a prototype denoted by p T u which is the prototype p u of the original map, u being the unit located at T -1 (k u 1 , k u 2 ). When comparing two maps, the mean of the squared distances (in R d ) between the prototypes of the two maps that are located at the same position is calculated. For two maps M and M , with respective prototypes (p u ) u and (p u ) u , we define a distance between two maps as the distance between their respective prototypes positionned at the same coordinates:

D (M, M ) = 1 m 2 m 2 u=1 p u -p u 2 . ( 1
)
The best transformation between the current fused map and the next map to be fused is chosen according to this distance. The two maps are then fused using the optimal transformation before they are merged, as described in Algorithm 1.

The optimal transformation is found by computing the distance between the maps to be fused, T * b (M b ), and a reference map, which can be the first of the Algorithm 1 Optimal transformation

Initialization M * ,1 ← M 1 for b : 2 → B do Optimal transformation T * b := arg min T ∈T D M * ,b-1 , T (M b ) Fusion between M * ,b-1 and T * b (M b ). Provides: M * ,b end for Return M * := M * ,B
list, M1 , for instance 1 . The fusion between the map is performed as suggested in [START_REF] Vrusias | Distributing SOM ensemble training using grid middleware[END_REF] by averaging the prototypes located at the same position:

∀ u = 1, . . . , m 2 , p * u := 1 B B b=1 p b,T u . (2)
In the method described in the previous section, all maps are fused in an arbitrary order. However, as pointed out in [START_REF] Pasa | Fusion of Kohonen maps ranked by cluster validity indexes[END_REF], the maps may have very different qualities and may also be very different: merging a very peculiar map with a poor quality might lead to deterioration of the the results instead of improving them. In this section, two strategies are presented to leverage this problem.

The first one uses a measure of quality of the maps and first rank the maps from the one with the best quality to the one with the worse quality: M (1) , . . . , M (B) . Standard quality measures for SOM can be used to perform this ranking [START_REF] Polzlbauer | Survey and comparison of quality measures for self-organizing maps[END_REF]: i) the quantization error (QE),

m 2 u=1 i: xi∈C * u x i -p * u 2
, which is a clustering quality measure, disregarding the map topology; ii) the topographic error (TE) which is the simplest of the topographic preservation measure: it counts the ratio of second best matching units that are in the direct neighborhood on the map of the best matching units for every (x i ) i . However, for small maps and relatively simple problems, this measure has a small variability and can lead to many equally ranked maps. Therefore, another approach is introduced to make a trade-off, while ranking the maps, between clustering and topographic qualities: the average rank of the maps is computed as:

r b = r b quanti + r b topo 2 (3)
where r b quanti is the rank of the map M b according to its quantization error (the best map is ranked first) and similarly for r b topo with the topographic error and the maps were finally ranked by increasing order of (r b ) b .

Taking advantage of this ordering of the maps, the previous method can be modified using two different strategies:

1. the similarity strategy: following an idea similar to [START_REF] Pasa | Fusion of Kohonen maps ranked by cluster validity indexes[END_REF], the maps are merged by similarity: the merging process is initialized with the best map: M * ,1 ← M (1) . Then, this map is merged only with the maps that resemble this reference map. To do so, a simple ascending hierarchical clustering is performed between the maps (T * b (M b )) b=1,...,B , with (T * b ) b obtained by comparison with the reference map M 1 . This clustering is based on the distance introduced in (1) and the hierarchical tree is cut using the method described in [START_REF] Langfelder | Defining clusters from a hierarchical cluster tree: the dynamic tree cut package for R[END_REF]. Finally, the maps in the same cluster as M (1) are fused to M * ,1 ; 2. the ordering strategy: an alternative approach is performed sequentially by merging the maps by increasing rank M (1) , M (2) , . . . The merging process is stopped at M (B ) with B ≤ B (and usually B < B) when the quality of the fused map M * ,B would not increase anymore by merging it with M (B +1) (actually, two strategies are investigated: stopping when the quality measure is not increasing or stopping when the quality measure has not increased for the last 5%B fused maps).

Simulations

Methodology. In all the simulations, B = 100 maps are generated using the standard SOM. The optimal B has not been investigated in this paper and the number of fused maps was simply taken large enough so that the fusion makes sense. All maps were built with approximately m = n 10 units and 5 × n iterations of the stochastic algorithm and equipped with a Gaussian neighborhood controlled with the Euclidean distance between units on the grid. The size of the neighborhood was progressively decreased during the training. All simulations have been performed using the R package SOMbrero2 . The 100 maps are then fused using one of the strategies described below and the performance of the different methods are finally assessed using various quality criteria for the resulting maps M * : i) two criteria already mentioned in Section 3 that are standard to measure the quality of the SOM: i) QE and TE; ii) a criterion which uses the ground truth, when available (i.e., an a priori group for the observations), the normalized mutual information (NMI) [START_REF] Danon | Comparing community structure identification[END_REF] between the unit of the map and the a priori group. This criterion quantifies the resemblance between the a priori group and the clustering provided by the SOM (it is comprised between 0 and 1, a value of 1 indicating a perfect matching between the two classifications). Note that this criterion must be interpreted with care because if the a priori groups are split between several units of the map, each of these units being composed of one group only (which is expected for SOM results), the criterion can be lower than when the groups are split between less units which are all composed of several groups (which would be a less expected result). Thus, this criterion has to be interpreted only together with the QE and the TE values.

The performance of the method is also assessed in term of stability. It is expected that several runs of one aggregating method give similar (thus stable) results. This stability is estimated in terms of: i) the distance between two final maps obtained from two different runs of the same method. If M * and M * are two maps, the quantity D(M * , T * ( M * )), where D is defined as in (1) and T * := arg min T ∈T D(M * , T ( M * )), is computed. This gives an estimation of the dissemblance between two maps from the prototype (hence the topological) perspective. If calculated over 250 different final maps, this quantity helps to quantify the stability of the final prototypes provided by a given aggregation method; ii) the NMI between the final classes of two final maps obtained from two different runs of the same method. This gives an estimation of the dissemblance from the clustering perspective for a given aggregation method.

250 fusions for each method are performed using the methodology described above. This permits to compute average quality as stability criteria as well as to have an overview of the distribution of these criteria when the method is repeated.

Compared methods. The comparisons performed in this section aim at comparing our approach to existing ones (which are described in Section 2) as well as to investigate several options of the method (as discussed in Section 3). First, our method, which merges several maps obtained from several initialization states, is compared to the standard bagging approach, in which several maps are trained from bootstrap samples from the similar initialization states. More precisely, bootstrap strategies are:

the method denoted by B-Rand, which uses a common random initialization to learn B = 100 maps from 100 bootstrap samples coming from the original data set. Then, the prototypes that are positioned at the same coordinates, are averaged to obtain the final map M * (as suggested in [START_REF] Vrusias | Distributing SOM ensemble training using grid middleware[END_REF]); the method denoted by B-PCA, which uses a common PCA initialization to learn B = 100 maps from 100 bootstrap samples coming from the original data set (as suggested by [START_REF] Petrakieva | Bagging and bumping self organising maps[END_REF]). The PCA initialization consists of initializing the prototypes by regularly positioning them along the coordinates of the projection of the data set on the first two axis of the PCA. Then, the prototypes that are positioned at the same coordinates, are averaged to obtain the final map M * ; the method denoted by B-Seq, which uses a sequential initialization of the B = 100 maps: the first map is initialized randomly and trained with a bootstrap sample and the b-th map is initialized with the final prototypes of the (b -1)-th map and trained with another bootstrap sample. Finally, the final map M * , is obtained by averaging the prototypes of the B = 100 maps, that are positioned at the same coordinates, as suggested in [START_REF] Baruque | Fusion methods for unsupervised learning ensembles[END_REF].

These strategies are compared with our method and its bootstrap version, respectively denoted by RoSyF (for "Rotation and Symmetry Fusion") and B-RoSyF. RoSyF learns B = 100 maps, each from a different random initial state and using the whole data set (x i ) i=1,...,n and B-RoSyF learns B = 100 maps from 100 bootstrap samples coming from the original data set. Finally, we also compare RoSyF with the approach consisting in selecting only one map from the B maps, the map supposed to be the best for instance. More precisely, using the B = 100 maps generated during the training of the RoSyF method, we selected one of the B = 100 maps i) randomly (this method is denoted by Best-R), ii) with the smallest QE (this method is denoted by Best-QE or iii) with the smallest TE (this method is denoted by Best-TE).

Datasets and results. This section compares the results obtained on two datasets coming from the UCI Machine Learning Repository3 as available in the R package mlbench 4 . More precisely, the data "Glass" (n = 214, d = 10 and 7 a priori groups) [START_REF] Towell | Interpretation of artificial neural networks: mapping knowledge-based neural networks into rules[END_REF] and the data "Vowel" (n = 990, d = 10 and 11 a priori groups) [START_REF] Niranjan | Neural networks and radial basis functions in classifying static speech patterns[END_REF] are used. The SOM parameters are set to m = 5 and 1 000 iterations for "Glass" and m = 10 with 5 000 iterations for "Vowel". The different strategies, and especially the relevance of using different initial states instead of different bootstrap samples with the same initialization, is evaluated. The results are provided in Table 1.

First, note that for almost all quality criteria and datasets, RoSyF obtain better results than the methods based on different bootstrap samples (all differences are significant according to Wilcoxon test, risk 5%). B-RoSyF slightly deteriorates RoSyF performances. [START_REF] Cottrell | A statistical tool to assess the reliability of self-organizing maps[END_REF][START_REF] De Bodt | Statistical tools to assess the reliability of self-organizing maps[END_REF] reported that the SOM algorithm is highly insensitive to initialization if run on the same data set as compared to what is obtained if bootstrap samples are used. However, it seems that the quality of the aggregated map is much better when different initial states are used on the same data set rather than different bootstrap samples with a common initial state, whatever this initial state is. Second, the TE obtained by RoSyF is always the lowest, just after the one obtained by Best-TE (which always selects the map with the lowest TE) but with a better QE and a better NMI. Again, all these differences are significant according to Wilcoxon tests (risk: 5%). On a clustering quality point of view, RoSyF is the method that obtains the second lowest quantization error, just after Best-QE which is designed to select the map with the lowest QE. Also, from a classification point of view, its performance is also very good: in average, RoSyF ranks first for the NMI criterion. Also note that all quality criteria have a low variability: the standard deviations is almost always the lowest: RoSyF is the method which has the best coefficient of variation (mean divided by the standard deviation) for all quality criteria. Table 2 (and Figure 1 for the dataset "Vowel") provides a comparison of the stability criteria. For this data set, RoSyF has the best stability, either in term of prototype stability (even though B-PCA and B-Seq also have a good prototype stability) and even more in term of class stability. These differences are significant according to Wilcoxon tests (risk: 5%). The results indicate that the method is indeed appropriate to improve the quality of the final map but also that it is very stable and gives very similar results if used several times, with different initializations of the prototypes and different training of the merged maps.

The relevance of stopping the merging process before all the maps have been fused has also been evaluated 5 . This comparison shows that there is only a small benefit in stopping the merging process before all maps have been used: most strategies lead to an highly deteriorated TE. Only stopping the training process when TE increases (TE-Inc) or based on the similarity strategy described in 3 are valid approaches in terms of quality criteria. However, a stability analysis shows that all these strategies strongly deteriorate the stability of the final map: merging all maps is the approach that provides the best stability, either in term of prototype comparison than in term of class comparison, except for TE-Inc which provides a slightly more stable clustering but very different prototypes. All these strategies use only few maps (less than 10 maps in average), except again TE-Inc which uses 89.4 maps in average for the "Glass" dataset and is thus very close to the maximum number of available maps (100). Actually, additional simulations (not shown for the sake of paper length) merging more than 100 maps proved that the stability increases with the number of fused maps (up to a certain number which was for our dataset between 500 to 1000 maps). A trade-off has thus to be found between computational time required to generate a large number of maps and stability of the results. This question is still under study.

Conclusion

Although most work on SOM ensembles are based on bootstrapping techniques, this paper presents an approach allowing to explore different initial states for the map. The method improves the stability of the fused map, both in term of prototypes and in terms of clustering. We are currently investigating how to choose an optimal number B of maps to fuse as well as weighting schemes based on various quality criteria: this approach is already promising to improve the results, especially the stability of the final map.

Fig. 1 .

 1 Fig. 1. Normalized mutual information (NMI) between pairs of clusterings obtained from the 250 final maps generated by the different approaches.

Table 1 .

 1 Method performance comparison (mean and standard deviation of different quality criteria; QE has been multiplied by 100)

		B-RandB-PCA B-Seq B-RoSyF RoSyF Best-R Best-QEBest-TE
					"Glass"
	mean QE	855.10 855.93 854.97 609.84	597.81 595.09 560.69 593.68
	sd QE	10.30	9.43	9.24	23.10	9.82	15.52	5.45	13.96
	mean TE	11.95% 12.42% 11.77% 0.01%	0.01% 0.10% 0.04%	0.00%
	sd TE	6.09% 6.53%	6.45% 0.04%	0.07% 0.24% 0.17%	0.00%
	mean NMI 15.80% 15.77% 16.00% 18.92%	17.86%15.64% 16.37% 15.87%
	sd NMI	3.38% 3.15%	3.30% 2.09%	1.38% 2.20% 2.03%	2.21%
					"Vowel"
	mean QE	847.57 847.73 847.91 550.78	545.88 547.44 531.30 548.23
	sd QE	11.82	10.88	11.63	5.18	1.01	7.10	2.39	6.72
	mean TE	5.89% 6.06%	5.80% 0.07%	0.07% 0.19% 0.20%	0.00%
	sd TE	3.62% 3.46%	3.37% 0.10%	0.08% 0.14% 0.14%	0.00%
	mean NMI 7.11% 6.76%	7.03% 9.47%	9.57% 9.64% 9.53%	9.53%
	sd NMI	1.44% 1.37%	1.49% 0.12%	0.11% 0.66% 0.54%	0.72%

Table 2 .

 2 Method stability comparison (mean and standard deviation of different stability criteria; D has been multiplied by 10 000)

		B-RandB-PCA B-Seq B-RoSyF RoSyF Best-R Best-QEBest-TE
					"Glass"
	mean D	70.85	67.22	67.06 149.65	67.07 2047.14 1302.27 1581.49
	sd D	38.62	32.32	31.24 335.14	310.74 1557.08 1170.39 1186.28
	mean NMI 64.77% 65.60% 65.88% 83.54%	87.47%49.15% 54.41% 49.86%
	sd NMI	6.37% 6.32%	6.23% 5.83%	5.11% 10.81% 9.57% 10.26%
					"Vowel"
	mean D	59.89	61.33	59.21	15.30	11.07 681.87 535.32 716.81
	sd D	31.19	33.33	31.42	5.77	3.87	275.23 185.06 343.41
	mean NMI 57.32% 56.83% 57.70% 90.83%	92.39%72.53% 74.94% 72.11%
	sd NMI	5.32% 5.21%	5.20% 1.59%	1.33% 3.29%	2.66%	3.37%

The current fused map, M * ,b-1 has also been used as a reference map, with no difference in the final result. Using M 1 is thus a better strategy, because optimal transformation can be computed in parallel.

http://cran.r-project.org/web/packages/sombrero, version 1.0.

http://archive.ics.uci.edu/ml

http://cran.r-project.org/web/packages/mlbench

For the sake of paper length, detailed results are not reported but only described.


Steve Della Mora

Lapo Boschi
email: lapo.boschi@upmc.fr

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

 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 leastsquares 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 viceversa 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., [START_REF] Becker | A comparison of tomographic and geodynamic mantle models[END_REF][START_REF] Bull | Global scale models of the mantle flow field predicted by synthetic tomography models[END_REF]. 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., [START_REF] Doornbos | Regions of seismic wave scattering in the earth's mantle and precursors to PKP[END_REF][START_REF] Haddon | Evidence for scattering of seismic pkp waves near the mantle-core boundary[END_REF][START_REF] Cormier | Anisotropy of heterogeneity scale lengths in the lower mantle from PKIKP precursors[END_REF][START_REF] Margerin | Multiple scattering of high-frequency seismic waves in the deep earth: Pkp precursor analysis and inversion for mantle granularity[END_REF][START_REF] Garcia | Statistical study of seismic heterogeneities at the base of the mantle from pkp differential traveltimes[END_REF].

With this study we explore a "stochastic" approach alternative to tomography, introduced by [START_REF] Gudmundsson | Stochastic analysis of global traveltime data: mantle heterogeneity and random errors in the isc data[END_REF] (hereafter GDC90) and [START_REF] Davies | Spectra of mantle shear wave velocity structure[END_REF] 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 spherical-harmonic 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-of-magnitude 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. [START_REF] Bunge | Effect of depth-dependent viscosity on the planform of mantle convection[END_REF][START_REF] Mégnin | Imaging 3-d spherical convection models: what can seismic tomography tell us about mantle dynamics?[END_REF][START_REF] Yoshida | Mantle convection with longest-wavelength thermal heterogeneity in a 3-d spherical model: Degree one or two?[END_REF][START_REF] Van Heck | Planforms of self-consistently generated plates in 3d spherical geometry[END_REF][START_REF] Foley | Generation of plate-like behavior and mantle heterogeneity from a spherical, viscoplastic convection model[END_REF][START_REF] Dziewonski | Mantle anchor structure: An argument for bottom up tectonics[END_REF][START_REF] Nakagawa | The influence of morb and harzburgite composition on thermo-chemical mantle convection in a 3-d spherical shell with self-consistently calculated mineral physics[END_REF][START_REF] Davies | Reconciling dynamic and seismic models of Earth's lower mantle: the dominant role of thermal heterogeneity[END_REF]. 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 travel-time 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.

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) δp(r, θ, ϕ) (with r, θ, ϕ radius, colatitude and longitude, respectively) and anomalies δt in seismic travel-time.

Neglecting non-linear effects, this relationship has the general form

δt = V K(r, θ, ϕ)δp(r, θ, ϕ)dV, (1)
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 δt. 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 non-zero on the ray path (traced in the reference model), and zero everywhere else. If some form of finite-frequency theory is used, K becomes more complicated (e.g., [START_REF] Peter | Surface wave tomography: global membrane waves and adjoint methods[END_REF][START_REF] Peter | Tomographic resolution of ray and finitefrequency methods: a membrane-wave investigation[END_REF]. When the assumption of linearity is dropped, e.g. if we care about multiple-scattering, then no function K can be defined, and eq. (1) ceases to be valid. Typically, eq. ( 1) is used to set up an inverse problem with δp(r, θ, ϕ) as the unknown, and a set of observations of δt as data.

Eq. ( 1) can be re-written for each observed value of δt, all of them with their corresponding kernel function K. δp is then expressed as a sum of unknown coefficients multiplied by some known "basis functions", e.g.

δp(r, θ, ϕ) = L l=0 l m=-l N n=1 A lmn Y lm (θ, ϕ)R n (r), (2)
with Y lm denoting the real scalar spherical harmonic of degree l and order m (e.g., [START_REF] Dahlen | Theoretical Global Seismology[END_REF] 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 mixeddetermined inverse problem with unknown coefficients A lmn , while all other quantities in (1) can be calculated.

The main idea of GDC90 and [START_REF] Davies | Spectra of mantle shear wave velocity structure[END_REF] is to set up an inverse problem whose unknowns are not the coefficients A lmn , but the spectral power per unit area

Q ln = 1 2l + 1 l m=-l A 2 lmn . (3)
This is achieved by first defining an approximately equal-area grid spanning the Earth's surface. All δt observations associated with sources/receivers lying in the same pair of grid cells ("bins") are grouped in a "summary ray" (e.g., [START_REF] Morelli | Topography of the core-mantle boundary and lateral homogeneity of the liquid core[END_REF] 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 Θ, range of epicentral distance (distance "bin", identified by its mean value ∆) and range of source depth (depth bin, identified by its mean value Z), a value of the variance of δt is calculated,

σ 2 (Θ, ∆, Z) = n S k=1 n k -1 n S k=1 n k n k i=1 [δt i -mean k (δt)] 2 n k -1 = n S k=1 n k -1 n S k=1 n k σ2 k (4) (GDC90)
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 (δt) is the mean of all measurements of δt within that bundle. The factor n k in (4) is introduced so that σ 2 is more strongly affected by summary rays formed by larger numbers of δt observations. For a given binning scheme (i.e., given values of Θ and Z), σ 2 is a function of the epicentral distance ∆. Through (4), the numerical values of σ 2 can be determined directly from the observations δt. 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 (Θ, ∆, Z) with n S > 10. For body waves, we only take into account rays with travel-time |δt| < 4s and ∆ < 100 • , bundles with n k > 4 and bins with n S > 4, in analogy with GDC90.

σ 2 as defined by eq. ( 4) can be thought of as the average, calculated over all ray bundles in the same (Θ, ∆, Z) bin, of the variance σ2 k of δt calculated within each ray bundle. In practice, after introducing the operator

E C = 1 n k -1 n k i=1
(...) (sum extended over all rays within a bundle) and the expected value operator

E = n S k=1 n k -1 n S k=1
n k (...) (sum over all bundles in the same (Θ, ∆, Z) bin), eq. ( 4) takes the more compact form

σ 2 (Θ, ∆, Z) = E E C (δt -E C (δt)) 2 .
(5) 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 δp 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).

3 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 (Θ, ∆, Z)) as a different realization of the same experiment.

After some algebra, eq. ( 5) can be written

σ 2 (Θ, ∆, Z) = E E C (δt) 2 -E E C (δt) 2 . ( 6
)
The second term of the latter expression can be rewritten

E E C (δt) 2 = E 1 A A δt(x)dx 2 , ( 7
)
where A is the area of the grid cell of radius Θ and δt(x) is defined by eq. ( 1). Inverting the order of the operators E and E C , and applying Fubini's theorem (e.g., [START_REF] Thomas | Calculus and Analytic Geometry[END_REF], we find

E E C (δt) 2 = 1 A 2 A A E δt(x 1 )δt(x 2 ) dx 1 dx 2 . ( 8
)
Let us now consider the other term in eq. ( 6),

E E C (δt) 2 = 1 A A E δt(x 1 )δt(x 1 ) dx 1 . (9) Since 1 A A dx 2 = 1, E E C (δt) 2 = 1 A A E δt(x 1 )δt(x 1 ) dx 1 1 A A dx 2 , ( 10
)
and again based on Fubini's theorem,

E E C (δt) 2 = 1 A 2 A A E δt(x 1 )δt(x 1 ) dx 1 dx 2 . (11) If we define ρ = |x 2 -x 1 |, then δt(x 1 ) = lim ρ→0 δt(x 2 )
, and

E E C (δt) 2 = 1 A 2 A A lim ρ→0 E {δt(x 1 )δt(x 2 )} dx 1 dx 2 . ( 12
)
Substituting ( 8) and ( 12) into (6), we find

σ 2 (Θ, ∆, Z) = 1 A 2 A A lim ρ→0 E {δt(x 1 )δt(x 2 )} dx 1 dx 2 + - 1 A 2 A A E δt(x 1 )δt(x 2 ) dx 1 dx 2 (13)
In section 3.1 we shall show that, in the assumption of Gaussian, stationary and isotropic slowness perturbations δp, the expression E {δt(x 1 )δt(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 ρ = |x 1 -x 2 |, and not of x 1 and x 2 themselves. In section 3.2 we shall apply the ray-theory 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 eq. ( 13), leading to a relatively simple expression for σ 2 in terms of the harmonic r). Such equation constitutes the basis of GDC90's and our formulation of the inverse problem.

spectrum of the Earth, Q l (r) = N n=1 Q ln R n (
3.1 Relation between delay-time variance within ray bundles, and the harmonic spectrum of Earth's structure

We next show how E {δt(x 1 )δt(x 2 )} can be written in terms of the Earth's spectral coefficients Q l (r). Based on eq. ( 1),

E {δt(x 1 )δt(x 2 )} = E V V K 1 (r 1 )K 2 (r 2 )δp(r 1 )δp(r 2 )dr 1 dr 2 , (14)
with r i = (r 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

E {δt(x 1 )δt(x 2 )} = V V K 1 (r 1 )K 2 (r 2 )E {δp(r 1 )δp(r 2 )} dr 1 dr 2 . ( 15
)
Let us focus on the term E {δp(r 1 )δp(r 2 )} at the right-hand side of this expression. Using eq. (2), and denoting for simplicity

A lm (r) = N n=1 A lmn R n (r), E {δp(r 1 )δp(r 2 )} = l,m p,q E {A lm (r 1 )A pq (r 2 )} Y lm (θ 1 , ϕ 1 )Y pq (θ 2 , ϕ 2 ). ( 16
)
We make at this point the important assumption, consistent with GDC90, that the fluctuations of δp are described, at any radius r within the mantle, by a Gaussian stochastic process, or in other words that, if we apply a shift ∆r to the slowness perturbation map δp(r), the correlation between δp(r) and δp(r + ∆r) quickly drops to zero with growing ∆r.

This property of δp(r) implies that E {δp(r 1 )δp(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

E {δp(r 1 )δp(r 2 )} = f (r, |r 2 -r 1 |, ρ), (17)
with ρ the angular horizontal distance between r 1 and r 2 :

cos(ρ) = cos(θ 1 ) cos(θ 2 ) + sin(θ 1 ) sin(θ 2 ) cos(ϕ 1 -ϕ 2 ). ( 18
)
We next write f (r, |r 1 -r 2 |, ρ) as a sum of Legendre polynomials P l (cos ρ) (e.g., [START_REF] Dahlen | Theoretical Global Seismology[END_REF] with coefficients (2l + 1)c l (r, |r 2 -r 1 |)/4π, and

E {δp(r 1 )δp(r 2 )} = l c l (r, |r 2 -r 1 |) 2l + 1 4π P l (cos ρ). (19)
The factor (2l + 1)/4π allows to simplify eq. ( 19) after application of the addition theorem

2l + 1 4π P l (cos ρ) = m Y lm (θ 1 , ϕ 1 )Y lm (θ 2 , ϕ 2 ) ( 20
)
(e.g., Dahlen & Tromp, 1998, eq. B.74), resulting in the expression

E {δp(r 1 )δp(r 2 )} = l,m c l (r, |r 2 -r 1 |) Y lm (θ 1 , ϕ 1 )Y lm (θ 2 , ϕ 2 ). ( 21
)
Let us equate the alternative expressions for E {δp(r 1 )δp(r 2 )} found at the right-hand side of eqs. ( 16) and ( 21):

l,m c l (r, |r 2 -r 1 |) Y lm (θ 1 , ϕ 1 )Y lm (θ 2 , ϕ 2 ) = l,m p,q E {A lm (r 1 )A pq (r 2 )} Y lm (θ 1 , ϕ 1 )Y pq (θ 2 , ϕ 2 ). ( 22
)
It follows from ( 22) and the orthogonality of Y lm that

c l (r, |r 2 -r 1 |) Y lm (θ 2 , ϕ 2 ) = p,q E {A lm (r 1 )A pq (r 2 )} Y pq (θ 2 , ϕ 2 ), (23)
and from ( 23) and, again, the orthogonality of Y lm that

c l (r, |r 2 -r 1 |) = E {A lm (r 1 )A lm (r 2 )} . ( 24
)
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 eq. ( 24) takes the form

c l (r, |r 2 -r 1 |) = Q l (r)c(|r 2 -r 1 |). (25)
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 eq. ( 19), and the resulting expression into (15), to find

E{δt(x 1 )δt(x 2 )} = 1 4π V V K 1 (r 1 )K 2 (r 2 ) l [(2l + 1)c (|r 1 -r 2 |) Q l (r)P l (cos ρ)] dr 1 dr 2 , (26)
a direct, linear relation between variance of δt within a ray bundle, and the Earth's harmonic spectrum Q l . 28). ρ is the horizontal distance between (θ 1 , φ 1 ) and (θ 2 , φ 2 ) introduced in eq. ( 26).

3.2 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 ray-theory 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, eq. ( 26) can be rewritten

E {δt(x 1 )δt(x 2 )} ≈ 1 4π ray 1 ray 2 l (2l + 1)Q l (r) c (|r 1 -r 2 |) P l (cos ρ)ds 1 ds 2 . ( 27
)
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, page 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 τ between r 1 and r 2 equals

τ = d 2 + s 2 2 ( 28
)
(see Figure 1). Then, for a function g(r 1 , r 2 ) that depends only on the distance τ between r 1 and r 2 ,

ray 1 ray 2 g(r 1 , r 2 )ds 1 ds 2 = 2 ray 1 ∞ d τ √ τ 2 -d 2 g(τ )ds 1 dτ. ( 29
)
The integral over τ in eq. ( 29) is particularly easy if g(τ ) = g 0 exp(-τ 2 /α 2 ) for some α,

implying ray 1 ray 2 g(r 1 , r 2 )ds 1 ds 2 = √ π ray 1 αg(d)ds 1 . (30)
Said x 1/2 the value of d such that g(x 1/2 ) = g 0 exp(-x 2 1/2 /α 2 ) = 1/2g 0 (i.e., x 1/2 is the "half-width" of the Gaussian g), it can be shown that x 1/2 = α √ ln 2, and

ray 1 ray 2 g(r 1 , r 2 )ds 1 ds 2 = π ln 2 ray 1 x 1/2 g(d)ds 1 . (31)
This result can be applied to the double integral at the right hand side of eq. ( 27), assuming that its argument, a function of the distance ρ between points on ray 1 and ray 2 , be close to Gaussian. Then

E {δt(x 1 )δt(x 2 )} ≈ 1 4 √ π ln 2 ray 1 x 1/2 (r) l [(2l + 1)Q l (r)P l (cos(ρ))] ds. ( 32
)
Notice how the function c (|r1 -r2|) disappears between eq. ( 26) and eq. ( 32), once the assumption on the correlation between δp on the whole planet is assumed to be Gaussian. As already mentioned by GDC90, eq. ( 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 ρ 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 eq. ( 13), i.e.

A 2

A A E {δt(x 1 )δt(x 2 )} dx 1 dx 2 . In sections 3.1 and 3.2 we have rewritten the integrand E {δt(x 1 )δt(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 σ 2 as datum, we reduce analytically the double surface integral

A A in (13)
to a single, one-dimensional integral.

Cartesian case

Consider the integral

I = A dx 1 A f (x 1 , x 2 ) dx 2 , (33)
with A a circular surface of radius Θ, 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 ρ 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 eq. ( 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

I = Θ 0 sds 2π 0 dχ A f (s, χ, x 2 ) dx 2 . ( 34
)
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 ρ (length) and ψ (angle), defined with respect to the location x 1 , or (s, χ). By definition, ρ 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 ρ and ψ. ρ ranges between 0 and s + Θ. For each ρ, the interval of values of ψ for which (ρ, ψ) falls within A must be determined (and integrated over). If s + ρ < Θ, that is ρ < Θ -s, such interval is (0, 2π). If ρ > Θ -s, the length φ of the arc of ψ to be integrated upon can be determined using the cosine rule,

Θ 2 = s 2 + ρ 2 -2sρ cos(φ/2) φ = 2 cos -1 s 2 + ρ 2 -Θ 2 2sρ , (35)
(see Figure 2 for a visual explanation of s, ρ and Θ).

Eq. ( 34) now becomes

I = Θ 0 sds 2π 0 dχ Θ-s 0 ρdρ 2π 0 f (s, χ, ρ, ψ) dψ + Θ+s Θ-s ρdρ φ 0 f (s, χ, ρ, ψ) dψ , ( 36
)
where we have made use of the fact that f depends only on ρ, 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, χ, ρ, ψ) = f (ρ), and

I = Θ 0 2πsds Θ-s 0 2πρf (ρ)dρ + Θ+s Θ-s 2ρ cos -1 s 2 + ρ 2 -Θ 2 2sρ f (ρ)dρ . ( 37
)
Figure 2: Sketch of variables s, ρ and Θ in eq. ( 35). A is the surface where the integrals in eq. ( 33) are done, C is the centre of A and Θ and s are the same as in eq. ( 35).

To further simplify the expression for I, it is convenient to change the order of integration over s and ρ. The first term at the right hand side of ( 37)

I 1 = 4π 2 Θ 0 sds Θ-s 0 ρf (ρ)dρ = 4π 2 Θ 0 ρf (ρ)dρ Θ-ρ 0 sds. ( 38
)
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 ρ between Θ -s and Θ, the other denoted I 3 , containing an integral over ρ between Θ and Θ + s. Changing the order of integration, we find

I 2 = 2π Θ 0 sds Θ Θ-s 2ρf (ρ) cos -1 s 2 + ρ 2 -Θ 2 2sρ dρ = 2π Θ 0 ρf (ρ)dρ Θ Θ-ρ 2s cos -1 s 2 + ρ 2 -Θ 2 2sρ ds, (39)
I 3 = 2π Θ 0 sds Θ+s Θ 2ρf (ρ) cos -1 s 2 + ρ 2 -Θ 2 2sρ dρ = 2π 2Θ Θ ρf (ρ)dρ Θ ρ-Θ 2s cos -1 s 2 + ρ 2 -Θ 2 2sρ ds. ( 40
)
If one now combines I = I 1 +I 2 +I 3 , eq. ( 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

x cos -1 x 2 + a bx dx = 1 2 x 2 cos -1 x 2 + a bx + - 2 √ -z + (b 2 -4a) tan -1 b 2 -2x 2 -2a 2 √ -z 8 ( 41
)
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,

I = 4πΘ 2 2Θ 0 cos -1 ρ 2Θ - ρ 2Θ 1 - ρ 2Θ 2 ρf (ρ)dρ. ( 42
)

Spherical case

To derive eq. ( 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, eq. ( 42) is replaced by eq. ( 27) of GDC90

I = A 2 2Θ 0 w(Θ, ρ)f (ρ)dρ, (43)
where

w(Θ, ρ) =        π -4 cos Θ cos -1 α + cos -1 β 1 + cos -1 β 2 2π(1 -cos Θ) 2 sin ρ if 0 < ρ < Θ π -4 cos Θ cos -1 α + sin -1 β 1 -sin -1 β 2 2π(1 -cos Θ) 2 sin ρ if Θ < ρ < 2Θ, (44)
and α, β 1 , β 2 are defined

α = cos Θ(1 -cos ρ) sin Θ sin ρ , (45)
β 1 = (1 -cos ρ) [1 + cos ρ -cos Θ(1 + cos Θ)] (1 -cos Θ) sin Θ sin ρ , (46)
and

β 2 = (1 -cos ρ) [1 + cos ρ + cos Θ(1 -cos Θ)] (1 + cos Θ) sin Θ sin ρ . ( 47
)
3.4 Relation between variance σ 2 , and the harmonic spectrum of Earth's structure

We have shown in section 3.

1 that E {δt(x 1 )δt(x 2 )} is a function of the distance ρ = |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 δt(x 1 )δt(x 2 ) = f (ρ), the double integral over surface is reduced to only one integral over distance by means of a weight function w(Θ, ρ), for both the cases of flat and spherical Earth. Replacing f (ρ) with E δt(x 1 )δt(x 2 ) in eq. ( 43),

1 A 2 A A E {δt(x 1 )δt(x 2 )} dx 1 dx 2 = 2Θ 0 w(Θ, ρ)E {δt(x 1 )δt(x 2 )} dρ. ( 48
)
After substituting the term E δt(x 1 )δt(x 2 ) with its explicit expression in (32), eq. ( 48) becomes

1 A 2 A A E {δt(x 1 )δt(x 2 )} dx 1 dx 2 = 1 4 √ π ln 2 2Θ 0 w(Θ, ρ) ray 1 x 1/2 (r) l [(2l + 1)Q l (r)P l (cos(ρ))] ds dρ. ( 49
)
Let us now consider the first term at the right-hand side of eq. ( 13). Using eq. ( 32) with ρ = 0, we find

E δt(x 1 )δt(x 1 ) = lim ρ→0 E {δt(x 1 )δt(x 2 )} = 1 4 √ π ln 2 ray 1 x 1/2 (r) L l=0 (2l +1)Q l (r) ds, ( 50
)
where we have used the fact that P l (1) = 1, independent of l. Substituting ( 49) and ( 50) into (13),

σ 2 (Θ, ∆, Z) = L l=0 1 4 √ π ln 2 2Θ 0 w(Θ, ρ) 1 -P l (cos ρ) dρ ray x 1/2 (r)(2l + 1)Q l (r)ds (51)
where σ 2 (Θ, ∆, Z) can be evaluated from the data, and at the right-hand side everything but Q l (r) is known. Eq. ( 51) constitutes the basis of the inverse problem solved by GDC90.

4 Formulation of the inverse problem in a 2-D description of surface-wave propagation

Projection to two dimensions

In a JWKB ray-theory description of surface-wave propagation (e.g., [START_REF] Ekström | Measurements and global models of surface wave propagation[END_REF],

δt(x, ω) = - 1 [v(ω)] 2 S K(θ, ϕ, ω)δv(θ, ϕ, ω)dΩ, (52)
where δv(θ, ϕ, ω) denotes lateral heterogeneities in surface-wave phase velocity v at angular frequency ω, and K is the corresponding sensitivity kernel. δv(θ, ϕ, ω) can naturally be rewritten as a linear combination of real spherical harmonics

δv(θ, ϕ) = l,m A lm Y lm (θ, ϕ). ( 53
)
with constant A lm (no r-dependence). Here and in the following we shall neglect the dependence of δv, v, K and δt on ω; in practice, we shall always consider different frequencies separately.

We assume, as in the 3-D case, that the variations of δv(θ, φ) 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 meaningless, be removed. In a spherical reference frame, eq. ( 15) can be rewritten

E {δt(x 1 )δt(x 2 )} = 1 v 4 S S K 1 (θ 1 , ϕ 1 )K 2 (θ 2 , ϕ 2 )E {δv(θ 1 , ϕ 1 )δv(θ 2 , ϕ 2 )} dΩ 1 dΩ 2 . ( 54
)
In analogy with section 4.2, we rewrite E {δv(θ 1 , ϕ 1 )δv(θ 2 , ϕ 2 )} by either replacing δv 1,2

with their harmonic expansions,

E {δv(θ 1 , ϕ 1 )δv(θ 2 , ϕ 2 )} = l,m p,q E {A lm A pq } Y lm (θ 1 , ϕ 1 )Y pq (θ 2 , ϕ 2 ) ( 55
)
(which is the 2-D counterpart of eq. ( 16)), or by invoking the statistical properties of δv (Gaussian, stationary, isotropic), which allow us to repeat steps ( 17) through (21) (after dropping the r-dependence of f ), resulting in

E {δv(x 1 )δv(x 2 )} = l,m c l Y lm (θ 1 , ϕ 1 )Y lm (θ 2 , ϕ 2 ), (56)
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 eq. ( 26),

E {δt(x 1 )δt(x 2 )} = 1 4πv 4 S S K 1 (θ 1 , ϕ 1 )K 2 (θ 2 , ϕ 2 ) l (2l + 1)Q l P l (cos ρ)dΩ 1 dΩ 2 , (57)
with ρ denoting the angular distance between x 1 and x 2 . 59).

In the ray-theory approximation,

E {δt(x 1 )δt(x 2 )} = 1 4πv 4 ray 1 ray 2 L l=0 (2l + 1)Q l P l (cos ρ)ds 1 ds 2 . (58)
If, in analogy with section 3.2, rays are treated as parallel, then a distance d between the two rays can be defined, and

ρ = d 2 + s 2 2 ( 59
)
(see Figure 3 for a visual explanation of ρ, s 2 and d); eq. ( 58) then becomes

E {δt(x 1 )δt(x 2 )} = 1 4πv 4 ray π d L l=0 (2l + 1)Q l ρP l (cos ρ) ρ 2 -d 2 dsdρ. ( 60
)
The ρ-integral in eq. ( 60) is evaluated between d and π, because we shall only consider firstorbit surface-wave observations, hence ρ < π. Since Q l does not depend on ρ, eq. ( 60) is rewritten

E {δt(x 1 )δt(x 2 )} = 1 4πv 4 ray π d L l=0 (2l + 1)Q l ρ ρ 2 -d 2 P l (cos ρ)dsdρ. ( 61
)
The ρ-dependent portion of eq. ( 61) forms a convergent integral that can be solved analytically (Appendix A). We can swap summation and integration in (61), and define

γ l (d) = π d ρ ρ 2 -d 2 P l (cos ρ)dρ. ( 62
)
We show in Appendix A how to calculate γ l (d). Substituting eq. ( 62) into (61) and keeping in mind that, in the parallel-ray approximation, ρ = d, we are left with the compact expression

E {δt(x 1 )δt(x 2 )} = 1 4πv 4 L l=1 (2l + 1)Q l ray γ l (ρ)ds. (63)
4.2 Relation between variance σ 2 (Θ, ∆) and the harmonic spectrum of surface-wave phase velocity heterogeneity

As noted in section 4.1, eq. ( 13) is valid, in the same form, in both the body-wave and surface-wave cases. Making use of equations ( 61)-( 62), the second term at the right-hand side of ( 13) can be rewritten

E E C (δt) 2 = 1 4πv 4 A 2 L l=0 (2l + 1)Q l A A ray γ l (ρ)ds. ( 64
)
According to eq. ( 12), the first term at the right-hand side of (13) coincides with the limit of the second as ρ → 0. Let us focus first on the integrand at the right-hand side of ( 13):

E δt(x 1 )δt(x 1 ) = lim ρ→0 E δt(x 1 )δt(x 2 ) = 1 4πv 4 L l=0 (2l + 1)Q l ray lim ρ→0 γ l (ρ)ds. ( 65
)
Details about the calculation of lim ρ→0 γ l (ρ) are shown in Appendix A. Substituting (107) into (65), and the resulting expression, together with (64), into (13), we are left with

σ 2 (Θ, ∆) = 1 4πv 4 A 2 L l=0 (2l + 1)Q l A A ray lim ρ→0 [γ l (ρ)] -γ l (ρ)ds dx 1 dx 2 . ( 66
)
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(Θ, ρ), so that eq. ( 66) becomes

σ 2 (Θ, ∆) = 1 4πv 4 L l=0 (2l + 1)Q l 2Θ 0 w(Θ, ρ) ray lim ρ→0 [γ l (ρ)] -γ l (ρ)ds dρ. ( 67
)
Since γ l (ρ) does not depend on s,

σ 2 (Θ, ∆) = ∆ 4πv 4 L l=0 (2l + 1)Q l 2Θ 0 w(Θ, ρ) lim ρ→0 [γ l (ρ)] -γ l (ρ) dρ, ( 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.

5.1

The 2-D problem: surface-wave dispersion and phase-velocity spectrum

Discretization

We approximate the integral in eq. ( 68) with a discrete summation, to find

σ 2 (Θ, ∆) = ∆ 4πv 4 L l=0 (2l + 1)Q l M m=1 w(Θ, ρ m ) lim ρ→0 [γ l (ρ)] -γ l (ρ m ) δρ, (69)
with

δρ = 2Θ M ρ m = (m -1)δρ + δρ 2 = m - 1 2 2Θ M . (70)
We next discretize values of cell size Θ and angular distance ∆, and eq. ( 69) takes the form

σ 2 (Θ i , ∆ j ) = ∆ j 4πv 4 L l=0 (2l + 1)Q l M m=1 w(Θ i , ρ m,i ) lim ρ→0 [γ l (ρ)] -γ l (ρ m,i ) δρ i . ( 71
)
While the discretization of ∆ is straightforward (we simply subdivide the domain ∆ = 0 A one-to-one correspondence is established between couples (i, j) and the values of a single index n, and

• < ∆ < 180 • into 2 • intervals), that
σ 2 n = l,m ∆ n 4πv 4 w(Θ n , ρ m,n ) lim ρ→0 [γ l (ρ)] -γ l (ρ m,n ) δρ n (2l + 1)Q l , (73)
which, after denoting

D n = σ 2 n F nl = M m=1 ∆ n 4πv 4 (2l + 1)w(Θ n , ρ m,n ) lim ρ→0 [γ l (ρ)] -γ l (ρ m,n ) δρ n , (74)
takes the simpler form

D n = l l=0 F nl Q l (75)
or, in a tensorial notation,

D = F • Q (76)
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., [START_REF] Snedecor | Statistical Methods[END_REF]. 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 σ2 k as defined in (4), that is (e.g., [START_REF] Bevington | Data reduction and error analysis for the physical sciences[END_REF])

C nn = n S k=1 n k n S k=1 n k 2 - n S k=1 n 2 k n S k=1 n k σ2 k -σ 2 n 2 . ( 77
)
Eq. ( 76) then becomes

C -1 2 • D = C -1 2 • F • Q, ( 78
)
which can be solved in least-squares sense to find the spectrum Q l . The tensor

C 1 2 is defined so that C 1 2 • C 1 2 = C. Since C is diagonal, so are C 1 2 and C -1 2 , with diagonal entries C 1 2
nn and C

-1 2 nn , respectively.

Least-squares solution and norm damping

We systematically discretize (Θ, ∆) so that N > L, where N is the largest possible value of n. Eq. ( 78) is then an overdetermined problem which admits the least-squares solution

Q = C -1 2 • F T • C -1 2 • F -1 • C -1 2 • F T • C -1 2 D = F T • C -1 • F -1 • F T • C -1 D (79)
(e.g., [START_REF] Trefethen | Numerical Linear Algebra[END_REF], where the superscript T denotes a transpose matrix, and we have made use of the fact that C -1 2 is diagonal and so

C -1 = C -1 2 • C -1 2 .
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 eq. ( 79) is regularized (e.g., [START_REF] Menke | Geophysical Data Analysis: Discrete Inverse Theory[END_REF]. As a regularization constraint, we impose that the norm of the solution be minimum. The least-squares formula (79) becomes

Q = F T • C -1 • F + λ 2 I -1 • F T • C -1 • D, (80)
with λ a regularization or "damping" parameter to be selected.

We solve eq. ( 80) by means of Cholesky factorization of F T • C -1 • F (e.g., [START_REF] Press | Numerical Recipes in Fortran 77: The Art of Scientific Computing[END_REF] (from now on LS) and of the non-negative least-squares algorithm (from now on NNLS)

of [START_REF] Lawson | Solving Least Squares Problems[END_REF]. 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 [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the l-curve[END_REF] to select an adequate value of λ: after defining the solution norm

ν(λ) = L l=0 Q (λ) l 2 (81)
(where the superscript λ identifies the solution found with the corresponding value of the damping parameter), we divide it by the norm of the undamped solution ν(0), thus defining the normalized norm

ν(λ) = 1 ν(0) L l=0 Q (λ) l 2 . ( 82
)
For each λ, we also define the solution misfit

ζ(λ) = N n=1 F • Q (λ) n -D n 2 N n=1 D 2 n , (83)
and we can build the L-curve plotting the couples (ν(λ), ζ(λ)). Our preferred value of λ is the one corresponding to maximum curvature of the L-curve [START_REF] Hansen | The use of the l-curve in the regularization of discrete ill-posed problems[END_REF].

The 3-D problem: body-wave travel times and the depth-dependent spectrum of the mantle

We next discretize eq. ( 51) to solve the original 3-D problem of GDC90. We first transform the ray integral over s to one over radius and find

σ 2 (Θ, ∆, Z) = L l=0 1 4 √ π ln 2 2Θ 0 w(Θ, ρ) [1 -P l (cos ρ)] dρ R R bot ds dr x 1/2 (r)(2l + 1)Q l (r) dr, ( 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 (∆,Z) (index i)

and Θ (index j), it follows from eq. ( 84) that

σ 2 (Θ j , (∆, Z) i ) = K i k=1 ds i dr k δr k L l=0 x 1/2 (r k )(2l + 1)Q l (r k ) 1 4 √ π ln 2 2Θ j 0 w(Θ j , ρ)[1 -P l (cos ρ)]dρ , (85)
where K i is the number of layers crossed by the ray associated to the i th bin and ds i dr k δr k is 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 eq. ( 85) reads

D = M • X, (86)
with

D p = σ 2 (Θ p , ∆ p , Z p ) M pq = ds p dr kq δr kq 1 4 √ π ln 2 2Θp 0 w(Θ p , ρ)[1 -P lq (cos ρ)]dρ (87)
X q = x 1/2 (r kq) (2l q + 1)Q lq r kq .

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 eq. ( 86) in a norm-damped least-squares sense, i.e.

X = M T • M + λ 2 I -1 • M T D, (88)
which we implement via NNLS, choosing λ according to the L-curve criterion. As opposed to the treatment of section 5.1.2 above, 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 σ 2 based on a set of surface-wave phase delays or body-wave travel-times, equations ( 80) and ( 88), respectively, provide the corresponding least-squares 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 surface-wave tomography maps are derived with the algorithm of [START_REF] Boschi | High-and low-resolution images of the earth's mantle: Implications of different approaches to tomographic modeling[END_REF]. 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 Θ-∆-Z bin is sampled by the data.

Surface-wave phase-velocity spectrum inversions

We apply our spectral inversion method to the fundamental-mode surface-wave dispersion database of [START_REF] Ekström | Measurements and global models of surface wave propagation[END_REF], focusing for brevity on phase delays of Rayleigh waves at a period of 50s (∼ 65, 000 observations) and of 100s Love waves (∼ 37, 000 observations).

Figures 4a and 5a

show the tomographic phase-velocity maps which we obtained from the same data, using a least-squares, ray-theory algorithm (e.g., [START_REF] Boschi | High-and low-resolution images of the earth's mantle: Implications of different approaches to tomographic modeling[END_REF]), 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 ∼1000km, well above that of the longest-period waves considered here (100s surface waves traveling at ∼4 km/s, hence wavelength ∼400km) and thus the physical resolution limit of imaging. Figures 4a and5a are in very good agreement with earlier results obtained from the same database, and so are the corresponding power spectra shown in Figures 4b and5b; see in particular the maxima at degrees 2 and 5, corresponding to the ocean-continent signature, which, in this period range, e.g. [START_REF] Carannante | Databases of surface wave dispersion[END_REF] have found to be a robust feature, independent of the technique used to measure surface-wave phase dispersion.

In Figures 4b and5b 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 eq. ( 88) to determine the depth-dependent spectrum of mantle P -wave velocity from the data set of direct P -wave travel times of [START_REF] Antolik | Global event location with full and sparse data sets using three-dimensional models of mantle P-wave velocity[END_REF] [START_REF] Boschi | High-and low-resolution images of the earth's mantle: Implications of different approaches to tomographic modeling[END_REF]. The "datum" σ 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 σ 2 are weighted according to the corresponding value of Θ 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 [START_REF] Antolik | J362d28: a new joint model of compressional and shear velocity in the earth's mantle[END_REF] provides summary rays with 2 • bins, no smaller value of Θ p are used in the calculation of σ 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 [START_REF] Boschi | High-and low-resolution images of the earth's mantle: Implications of different approaches to tomographic modeling[END_REF]; 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 section 6.1: at all mantle depths except for a few hundred km below the transition zone, the tomography spectrum has a clear, well known (e.g., [START_REF] Becker | A comparison of tomographic and geodynamic mantle models[END_REF][START_REF] Dziewonski | Mantle anchor structure: An argument for bottom up tectonics[END_REF]) maximum at degree 2; the "stochastic" spectrum is 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, ζ = 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. 89), as a function of harmonic degree l and bin size Θ.

7.1 Surface-wave phase velocity spectrum resolution

7.1.1 Sensitivity of σ 2 to Q l
Eq. ( 74) describes the relationship between the unknown spectral coefficients Q l , and the value of σ 2 associated with a (Θ, ∆) bin. On the basis of ( 74), let us introduce a sensitivity function

K nl = M m=1 (2l + 1)w(Θ n , ρ m,n ) lim ρ→0 [γ l (ρ)] -γ l (ρ m,n ) δρ n , (89)
so that the actual kernel relating σ 2 at Θ n with the spectral coefficient Q l coincides with the product ∆ n K nl .

We show in Figure 7 sensitivity K nl as a function of the size Θ n of geographic bins, and the harmonic degree l of the unknown spectral coefficients to be inverted for. Because in eq. ( 74) the epicentral distance ∆ n acts as a simple scaling factor, K nl alone fully describes the sensitivity of σ 2 to Q l . As a general rule, we see from Figure 7 that different harmonic degrees are constrained by values of σ 2 associated to systematically different geographic binning, i.e. large Θ n are needed to constrain the low-degree spectrum, while smaller Θ n serves to determine the higher-degree portion of the spectrum. It is immediately evident from Figure 7 that the averaged variance σ 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 σ 2 at degrees 3 and 90) for the LS solution for a dataset of Rayleigh waves at 50 s from [START_REF] Ekström | Measurements and global models of surface wave propagation[END_REF].

higher. At relatively high harmonic degrees, values of σ 2 associated to a broad range of Θ n values are sensitive 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 eq. ( 80) is

R = F T • C -1 • F + λ 2 I -1 • F T • C -1 • F (90)
(e.g., [START_REF] Menke | Geophysical Data Analysis: Discrete Inverse Theory[END_REF]. R is a measure of how well each model parameter (i.e. harmonic degree) is resolved in our inversions: the closer R is to the identity matrix, the higher the resolution. 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., [START_REF] Boschi | Measures of resolution in global body wave tomography[END_REF]. R depends on data coverage, but not on the data themselves. The value of the damping parameter λ, 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 50s fundamentalmode Rayleigh-wave data set of [START_REF] Ekström | Measurements and global models of surface wave propagation[END_REF], inverted in section 6.1. We implement eq. ( 90) via Cholesky factorization of F T •C -1 •F, using the value of λ selected according to the L-curve criterion when applied to the inversion of real, 50s 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 δt in the raytheory approximation, i.e., we implement eq. ( 52) integrating along the shortest great circle connecting source and receiver (no ray tracing). We then substitute the resulting synthetics into eq. ( 4) to find the synthetic σ 2 (Θ, ∆), 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 σ 2 (Θ, ∆) 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 σ 2 (Θ, ∆)).

Our first input model is constrained to have a simple monochromatic spectrum Q l = δ 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 50s 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 σ 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: loss-of-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 [START_REF] Boschi | High-and low-resolution images of the earth's mantle: Implications of different approaches to tomographic modeling[END_REF],

Data coverage and spectral resolution

To evaluate whether the lack of resolution anticipated in sections 7.1.1 and 7.1.2, and confirmed by the synthetic test of section 7.1.3, can be explained as an effect of inadequate (poor/nonuniform) data coverage, 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. [START_REF] Trampert | High resolution global phase velocity distributions[END_REF]. 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 σ 2 (Θ, ∆) 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.3 the misfit is high (ζ(λ) = 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.3 and 7.1.3. 7.2 Body-wave velocity spectrum resolution 7.2.1 Sensitivity of σ 2 to Q l Eq. ( 29) of GDC90 describes the sensitivity K 3D nl of σ 2 , for a given grid-size Θ, to the spectral coefficient Q l of mantle heterogeneity, i.e.

K 3D nl = 1 4 √ π ln 2 2Θn 0 w(Θ n , ρ) [1 -P l (cos ρ)] dρ. ( 91
)
The kernel K 3D nl is the 3-D counterpart of K nl as defined by eq. ( 89) above. Discretizing as explained in section 5.1.1,

K 3D nl = 1 4 √ π ln 2 M m=1 w(Θ n , ρ m,n ) [1 -P l (cos ρ m,n )] δρ n . ( 92
)
We plot K 3D nl in Figure 11 as a function of bin size Θ 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 σ 2 associated with large Θ n are sensitive to structure at relatively large harmonic degrees (e.g. for Θ n = 30 • , K 3D nl at l = 5 is about as large as K 3D nl at l = 40.

Because sensitivity of σ 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 σ 2 at large Θ to high degree Q l , is in good qualitative agreement with Figure 10b of GDC90.

Resolution matrix

The resolution matrix associated with eq. ( 88) is

R = M T • M + λ 2 I -1 • M T • M. (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 λ, 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 [START_REF] Antolik | J362d28: a new joint model of compressional and shear velocity in the earth's mantle[END_REF] (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 σ 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 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 [START_REF] Antolik | J362d28: a new joint model of compressional and shear velocity in the earth's mantle[END_REF] and adopting a linear relation between model anomalies and data (e.g., [START_REF] Boschi | High-and low-resolution images of the earth's mantle: Implications of different approaches to tomographic modeling[END_REF]. We calculate σ 2 (Θ, ∆, Z) via eq. ( 4), with M now defined by eq. ( 87). Consistently with section 7.1.3, we verify that this is approximately equivalent to substituting Q in eq. ( 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 ζ(λ) = 0.267 is good, which, together with the strong discrepancy between input and output spectra, indicates that the sensitivity of σ 2 (Θ, ∆, Z) to the Earth's spectrum is severely limited.

Discussion and conclusions

With this study we attempted to devise an 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 δp 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., [START_REF] Somerville | Cosmic variance in the great observatories origins deep survey[END_REF].

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., [START_REF] Carannante | Databases of surface wave dispersion[END_REF], who found highly correlated results from completely independent databases; for relatively short-period (∼ 30s) 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 [START_REF] Becker | A comparison of tomographic and geodynamic mantle models[END_REF], 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 [START_REF] Becker | A comparison of tomographic and geodynamic mantle models[END_REF] 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.4 we illustrate the results of a synthetic test involving 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 and12: 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 [START_REF] Mora | Beyond tomography: New insights in global seismic imaging[END_REF]; the nonlinear approach also allowed Della [START_REF] Mora | Beyond tomography: New insights in global seismic imaging[END_REF] to drop some of the approximations required here; yet, resolution remained equally poor (Della Mora, 2012, sec. 2.8). The very low sensitivity of σ 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.

acknowledgments

This work would not have been possible without the kind support and encouragement of Domenico Giardini and Martin Mai, who both also provided useful advice. Thorsten Becker helped with many insightful suggestions. We thank Ólafur Guðmundsson, Yanick Ricard, Jean Virieux and several anonymous reviewers for their detailed and insightful comments on the previous versions of this manuscript. [START_REF] Arfken | Mathematical Methods for Physicists[END_REF], we find

(2k -1)!! = 2 k Γ(k + 1 2 ) √ π (2k)!! = 2 k Γ(1 + k) (100) 2(l -k) -1 !! = 2 (l-k) Γ(l -k + 1 2 ) √ π 2(l -k) !! = 2 (l-k) Γ(l -k + 1),
where Γ(z) denotes the Euler's gamma function [START_REF] Whittaker | A course of modern Analysis[END_REF] Γ(z) = +∞ 0 e -t t z-1 dt. ( 101)

Replacing ( 100) into (98), the latter equation reduces to

P l (cos ρ) = 1 π l k=0 Γ(k + 1 2 ) Γ(k + 1) Γ(l -k + 1 2 ) Γ(l -k + 1) cos[(l -2k)ρ]. (102)
After deriving eq. ( 102) with respect to ρ and substituting the result into (95) (104)

Combining equations ( 94) and ( 103) we find the following expression for the integral (62): 112); right panel, the difference between γ l from the analytical evaluation of eq. ( 112) and the numerical one of eq. ( 108).

γ l (d) =        √ π 2 -d 2 if l = 0 (-1) l √ π 2 -d 2 + 1
where

δρ = π -d J ρ j = d + (j -1) + δρ 2 = d + δρ j - 1 2 = d + π -d J j - 1 2 , ( 111
)
so that the resulting expression of γ l (d) is

γ l (d) =                √ π 2 -d 2 if l = 0 (-1) l √ π 2 -d 2 + + 1 π l k=0 (l -2k) Γ(k+ 1 2 ) Γ(k+1) Γ(l-k+ 1 2 ) Γ(l-k+1) J j=1 ρ 2 j -d 2 sin[(l -2k)ρ j ]δρ if l = 0 . ( 112
)
The advantage of our approach in equation( 112) with respect to the integration of eq. ( 108) is that the term dP l (cos ρ) dρ oscillates much more intensively then sin [(l -2k) ρ] 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

Figure 1 :

 1 Figure1: Sketch of variables s 2 , τ and d in eq. (28). ρ is the horizontal distance between (θ 1 , φ 1 ) and (θ 2 , φ 2 ) introduced in eq. (26).

Figure 3 :

 3 Figure 3: Sketch of variables s 2 , ρ and d in eq. (59).

  of Θ 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 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 Θ in the range 3 • ≤ Θ ≤ 45 • , ordered by increasing value. Values are closely spaced (≤ 1 • ) up to Θ = 15 • , after which only Θ = 22.5 • , 30 • , 45 • is possible. Values of ∆ are regularly discretized from ∆ = 1 • to ∆ = 179 • with a sampling of 2 • .

Figure 4 :

 4 Figure 4: (a) Tomography map of 50s 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 :

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

  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

Figure 6 :

 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 7 :

 7 Figure7: Sensitivity K nl defined by eq. (89), as a function of harmonic degree l and bin size Θ.

Figure 8 :

 8 Figure 8: Resolution matrix from eq. (90) for the LS solution for a dataset of Rayleigh waves at 50 s from[START_REF] Ekström | Measurements and global models of surface wave propagation[END_REF].

Figure 9 :

 9 Figure 9: (a) Monochromatic random input model with Q l = δ 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: (a) Tomography model of 150s 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 :

 11 Figure 11: Same as Figure 7, but with the 3-D sensitivity function of eq. (92).

Figure 12

 12 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 :

 13 Figure13: 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.

  d 2 sin[(l -2k)ρ]dρ.In the following, we shall compact the notation by definingh l-2k (d) def = π d ρ 2 -d 2 sin[(l -2k)ρ]dρ.

Figure 14 :

 14 Figure 14: Validation of our analytical expression (112) for γ l ; (left panel), values of γ l from eq. (112); right panel, the difference between γ l from the analytical evaluation of eq. (112) and the numerical one of eq. (108).

A Analytical expression of γ l (d) and numerical validation A.1 Analytical integration of γ l (d)

We have seen in section 4.1 that an important part of the algorithm is the calculation of γ l (d). If l = 0 then eq. ( 62) becomes

If l = 0, after integration by parts of eq. ( 62),

The integral in (95) can be calculated with the approximated formula

where

Eq. ( 96) is a good approximation of (95) if W is large enough. The oscillations of dP l (cos ρ) dρ increase with l, so W must also increase with l if we want eq. ( 96) to be a good approximation of (95).

From [START_REF] Whittaker | A course of modern Analysis[END_REF],

Using the formula

Using the expression (102) for P l (cos ρ), this reduces to

A.2 Numerical implementation -Validation of our analytical expression for γ 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 γ l according to its analytically derived expression (105).

We validate our analytical integration of eq. ( 62), carried out in section A, by comparing its result (103) with the values of γ l found from eq. ( 94) (l = 0) and by direct numerical integration of eq. ( 95). We do not integrate eq. ( 62) numerically because it is singular at ρ = d. Substituting eq. ( 96) in (95),

In our approach, we substitute in eq. ( 105) the definition of h l-2k (d) of eq. 104, obtaining (110)