Franck Latallerie

Christophe Zaroli
email: c.zaroli@unistra.fr

Sophie Lambotte

Alessia Maggi

Analysis of tomographic models using resolution and uncertainties: a surface wave example from the Pacific

Keywords: Seismic tomography, Inverse theory, Surface waves and free oscillations, Pacific Ocean

teaching and research institutions in France or abroad, or from public or private research centers.

paths; other regions may be poorly constrained. This irregular data coverage causes the inverse problem to be under-determined: within poorly covered parts of a region, the model can change widely without changing how it fits the data. Which model we choose is, in some sense, arbitrary (Scales & Snieder 1997;[START_REF] Nolet | A Breviary of Seismic Tomography: Imaging the Interior of the Earth and Sun[END_REF]. Some tomographic studies make use of irregular or adapted meshes to account for the heterogeneous distribution of the data (e.g. [START_REF] Sambridge | Seismic tomography with irregular meshes[END_REF]. Many tomographic studies use ad-hoc regularization constraints to reduce the nonuniqueness of the model solution; they often minimize the model's complexity by damping the inversion (e.g. [START_REF] Nolet | A Breviary of Seismic Tomography: Imaging the Interior of the Earth and Sun[END_REF], Chapter 14) -hence the name damped-least-squares inversion. Such inversions may underestimate seismic velocity anomalies where data coverage is poor; less intuitively, they can also overestimate velocity anomalies where data coverage is highly uneven. In either case, we call such models locally biased (in the sense of 'averaging bias' effects, as reported by [START_REF] Zaroli | Toward Seeing the Earth's Interior Through Unbiased Tomographic Lenses[END_REF].

Half-a-century ago, Backus and Gilbert proposed an original approach to invert linear geophysical data designed to produce complete uncertainty information, optimize local resolution, and eliminate bias [START_REF] Backus | Numerical Applications of a Formalism for Geophysical Inverse Problems[END_REF][START_REF] Backus | The Resolving Power of Gross Earth Data[END_REF][START_REF] Backus | Uniqueness in the inversion of inaccurate gross Earth data[END_REF]. Instead of searching for a model that fits the data with some smoothness criterion, the Backus-Gilbert method searches for a model that has optimal resolution given the geometry and uncertainties of the data. [START_REF] Pijpers | Faster formulations of the optimally localized averages method for helioseismic inversions[END_REF], 1993) reformulated the method, applied it to helio-seismology, and named it Subtractive Optimally Localized Averages (SOLA). More recently, [START_REF] Zaroli | Global seismic tomography using Backus-Gilbert inversion[END_REF] introduced and adapted the SOLA method to seismic tomography problems with large numbers of parameters; he then applied this new inversion scheme to obtain a body wave tomographic model of the lower mantle.

The SOLA Backus-Gilbert method still suffers from a form of non-uniqueness, as model resolution and model uncertainties trade-off against each other. Despite not being designed for that purpose, the method also leads to models that fit the data. Note that SOLA tomography can deal with heterogeneous spatial distribution of the data and usually produces models with no averaging bias. Computationally, the SOLA formulation of the Backus-Gilbert approach remains tractable even when dealing with large-scale tomographic problems [START_REF] Zaroli | Global seismic tomography using Backus-Gilbert inversion[END_REF][START_REF] Zaroli | Toward Seeing the Earth's Interior Through Unbiased Tomographic Lenses[END_REF]). Lastly, it provides the full resolution and model uncertainties, which are the model statistics necessary to draw well-informed conclusions from the tomographic model.

In this study, we use the SOLA inversion to produce tomographic images of the Pacific upper mantle from path-averaged shear-wave velocity profiles obtained from surface-waves. After discussing how resolution and uncertainty trade-off in our models, we illustrate how to use them to perform meaningful comparisons of tomographic images with a geodynamic prediction. Taking as a reference a simple half-space cooling model [START_REF] Parker | Thermal Model of Ocean Ridges[END_REF], we point out anomalous regions in our tomographic model and we argue about their statistical significance and whether they are resolved or not.

DATA AND METHOD

Most surface-wave tomography studies that produce 3D shear-wave velocity models proceed in a sequence of two steps: either they first make 2D phase-or group-velocity maps at various frequencies (linear tomographic inversion), then combine them to obtain depth-dependent shear-wave velocity profiles at each point in the map (non-linear inversion) -examples are [START_REF] Ekström | A global model of Love and Rayleigh surface wave dispersion and anisotropy, 25-250 s: Global dispersion model GDM52[END_REF] or [START_REF] Liu | Global Rayleigh wave phase-velocity maps from finite-frequency tomography[END_REF]; or they first invert for 1D shear-wave profiles as a function of depth along each Tomographic analysis with resolution and uncertainties 5 source-station path (non-linear inversion), then combine these path-averaged velocity profiles into maps at each depth (linear tomographic inversion) -examples are [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF] or [START_REF] Isse | Surface wave tomography for the Pacific Ocean incorporating seafloor seismic observations and plate thermal evolution[END_REF]. Each sequence contains both a linear and a non-linear inversion. Only the linear inversion can be performed using the SOLA method.

Path-averaged velocity and uncertainty

As we aim to discuss how SOLA can improve the way we interpret seismic tomography images, we have chosen to apply the method on a pre-existing data-set of path-averaged shear-wave velocity profiles obtained from Rayleigh waves [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF]). These profiles range from 50 km to 450 km depth, and were obtained by inverting over 56 000 multi-mode surface-waves whose paths crossed the Pacific Ocean, using the Debayle (1999) automation of the [START_REF] Cara | Waveform inversion using secondary observables[END_REF] secondary observables method. After obtaining a shear-wave velocity profile for each path, [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF] gathered similar paths (those whose end-points were within 200 km of each other) into ∼15 000 clusters, obtained the average velocity profile of each cluster (Fig. 1a), and estimated the uncertainty of the velocities using the standard-deviation in each cluster. These clustered pathaveraged shear-wave velocities are the input data for our SOLA inversion.

The robustness of the analysis later in this paper relies on the reliability of the data uncertainty estimates. The approach of [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF] to estimate the data uncertainty captures the stochastic uncertainty: it measures the variability between adjacent paths within one cluster due, for example, to random measurement errors or to random errors in the location of the earthquakes in the cluster. However, it does not account for errors that would affect identically all paths in a cluster; for example, neglecting off-great circle paths due to some lateral heterogeneity or mislocating all earthquakes in the same direction due to lateral heterogeneities or network configuration. [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF] decrease at larger depths because paths along which long-period and higher-mode data were unavailable were damped toward the background model, a smoothed version of PREM [START_REF] Dziewonski | Preliminary reference Earth model[END_REF]. As all the path-averaged models tend toward the same values when the sensitivity decreases, the stochastic uncertainties also decrease.

Data uncertainties from

Since the [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF] approach considers only stochastic uncertainties, it underesti-

χ 2 red ( m) = 1 N N i=1 [( M j=1 G ij m j ) -d i ] 2 σ 2 d i , (1)
where m j is the j-th model parameter, M the number of model parameters, d i the i-th datum, N the number of data, G ij the elements of the sensitivity matrix of size

N × M such that M j=1 G ij m j
is the i-th predicted data, and σ d i the i-th data uncertainty. The values of χ 2 red at all depths are much greater than one (see Table 1). Two reasons may explain these values. First [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF] may have chosen an overdamped solution to ensure that all features in their model were interpretable. Second, the data uncertainties may be underestimated because they do not account for systematic bias. To obtain more reliable data uncertainty estimates, we assume that the final model of [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF] is a model that explains the data so that the χ 2 red values should be one, then we rescaled the data uncertainties by the factor required to bring χ 2 red = 1, i.e. we multiplied each data uncertainty by (χ 2 red ) 1/2 at the corresponding depth (see Table 1).

Regionalization with SOLA

We divide the Pacific region into independent depth layers, each of them divided into cells of 2 • in latitude and longitude. Each model parameter m j (j = 1, ..., M ) represents the shear-wave velocity in cell j. Each datum d i (i = 1, ..., N ) represents the average shear-wave velocity along the i-th source-receiver path. The sensitivity matrix G linearly relates the data to the velocities as follows:

d i = M j=1 G ij m j + n i , (2)
where n i denotes a noise term (we assume uncorrelated zero mean Gaussian noise). The tomographic problem is ill-posed because some regions are under-sampled and contradictions may arise from the noise in the data. Therefore the linear sensitivity matrix G is not invertible and we have to seek a 'generalized' inverse matrix G † , such that the k-th parameter estimate m k can be written as a linear combination of the data:

m k = N i=1 G †(k) i d i = N i=1 M j=1 G †(k) i G ij m j + N i=1 G †(k) i n i . (3)
Here k refers to the k-th grid cell,

G †(k) = (G †(k) i
) i=1,..,N represents the k-th row of the G † matrix (k may vary from 1 to M) and R = G † G is called the resolution matrix. We define the k-th averaging (or resolving) kernel as:

A (k) = (A (k) j = R (k) j /S j ) j=1,...,M , (4)
where

R (k) j = N i=1 G †(k) i
G ij are the elements of the k-th row of the resolution matrix. The factor S j is, in our case, the cell surface area; it appears because the model parameterization has cells with different sizes.

If the problem were well-posed, the sensitivity matrix would be invertible, the resolution matrix would be the identity matrix, and the model estimate m would equal the true model m. In practice, the model estimates are local averages of the true Earth. The k-th row of the resolution matrix R (k) (closely related to the k-th averaging kernel through the term S j ) relates m k to the true

parameters ( m k ≈ M j=1 R (k) j m j = M j=1 S j A (k) j m j ).
In the Backus-Gilbert approach, we want each averaging kernel to be centered and peak-shaped around the location of its corresponding model parameter.

The model estimate m k also contains the term G †(k) n that reflects propagation of data noise into model space. Note that smaller averaging kernel sizes lead to higher model uncertainties [START_REF] Menke | Geophysical Data Analysis: Discrete Inverse Theory[END_REF]). This can be understood intuitively: fewer data can be averaged within smaller averaging kernels, leading to fewer chances for data errors to cancel out. There is, therefore, a trade-off between good resolution and low model uncertainties.

The model uncertainty is defined statistically as:

σ m k = N i=1 G †(k) i σ 2 d i . (5)
In this study we assume uncorrelated zero-mean Gaussian data uncertainty. Thus the data covariance matrix is diagonal:

C d = diag(σ 2 d i ), i = 1, .
.., N . Note that it would be possible to consider non-diagonal data covariance matrix in the SOLA approach [START_REF] Pijpers | Faster formulations of the optimally localized averages method for helioseismic inversions[END_REF]) but we ignore off-diagonal terms for simplicity.

The key idea of the SOLA method is to specify an a priori target form,

T (k) = (T (k) j ) j=1,...,M ,
for each averaging kernel, A (k) . We specify some a priori information on the model resolution, which is fundamentally different from specifying a priori information on the model itself (as in data-fitting methods such as Bayesian methods and damped-least-squares). Rather than minimizing the spread (spatial extent) of each averaging kernel, A (k) , as in the original Backus-Gilbert approach, the SOLA variant aims to minimize the integrated squared difference between each averaging kernel, A (k) , and its corresponding target kernel, T (k) . A minimization problem directly computes each generalized inverse vector for each k-th parameter:

       arg min G †(k) ∈R N M j=1 S j (A (k) j -T (k) j ) 2 + η 2 σ 2 m k s.t. M j=1 S j A (k) j = 1, (6)
where η is a trade-off parameter to give more weight to the resolution or to the model uncertainties in the optimization. Equation 6is independent for each model parameter therefore efficient parallel computations are possible. We could adapt the trade-off parameter η for each model parameter but this would reduce the optimization of calculation; here we adapt only the target resolution. The constraint that the averaging kernels should be unimodular means that SOLA tomographic models should be free of averaging bias; a condition that cannot usually be guaranteed with inversion methods incorporating regularization [START_REF] Zaroli | Toward Seeing the Earth's Interior Through Unbiased Tomographic Lenses[END_REF].

The target resolution for a parameter is non-zero only in a circular region around that parameter's location; the circle's radius is computed from the path density:

r(ρ) = r max -(r max -r min ) log 10 (ρ) -log 10 (ρ min ) log 10 (ρ max ) -log 10 (ρ min ) , ( 7
)
where r is the target kernel radius, ρ is the path density, ρ min and ρ max are the minimal and maximal path densities, and r min and r max are the minimal and maximal target kernel radii. The logarithms lower the spatial variations of the target kernel radii compared to those of the path density. This smoothing is important to avoid tomographic images that would be difficult to interpret if nearby averaging kernels had widely different sizes. Other formulae could be chosen to reduce even further the spatial variability of the target kernel radii. The optimization problem (equation 6) involves the rows of the generalized inverse through the averaging kernels A (k) (equation 4) and model uncertainties σ m k (equation 5). It leads to a set of independent equations whose unknowns are the rows of the generalized inverse. Derivation of these equations from equation 6 can be found in [START_REF] Zaroli | Global seismic tomography using Backus-Gilbert inversion[END_REF]. These equations are solved using the LSQR algorithm of [START_REF] Paige | LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares[END_REF] as proposed by [START_REF] Nolet | Solving or resolving inadequate and noisy tomographic systems[END_REF] and [START_REF] Zaroli | Global seismic tomography using Backus-Gilbert inversion[END_REF].

Once the k-th row of the generalized inverse G †(k) has been computed, one obtains the associated model solution m k from equation 3, the model uncertainty σ m k from equation 5, and the averaging kernel A (k) from equation 4. For further details, the reader is referred to [START_REF] Zaroli | Global seismic tomography using Backus-Gilbert inversion[END_REF], [START_REF] Zaroli | Toward Seeing the Earth's Interior Through Unbiased Tomographic Lenses[END_REF][START_REF] Zaroli | Seismic tomography using parameter-free Backus-Gilbert inversion[END_REF].

RESULTS

In this section, we present the effects of varying the trade-off parameter η in equation 6, choose a value for η, describe in more detail the features of the resulting tomographic model, and illustrate how to take resolution and uncertainty into account when comparing a tomographic model to the theoretical predictions of a simple geodynamic model (half-space cooling). coverage. This is the well-known smearing effect, visible in the averaging kernel in the north-east of Hawaii, where paths have a predominant north-east to south-west orientation. In regions where the azimuthal path coverage is better, we do not see such artifacts. This is the case, for example, in the north-east of Japan where the averaging kernels are compact and nearly circular. Improving the resolution by lowering η makes the averaging kernels more circular and decreases the length and magnitude of the ramifications. The averaging kernels may have a small negative component, which is physically meaningless. A strong negative component would indicate that the target resolution was badly chosen and invalidate the results; in our case, it is small enough to be ignored.

A complete set of averaging kernels -one map for each model parameter -would fully represent the resolution of a SOLA inversion, but would be cumbersome and difficult to use when interpreting tomographic models. To simplify this information, we compute for each model parameter k a resolution length L k corresponding to the radius of a circle that contains 68% of the averaging kernel -a proxy to the standard deviation for a 2D Gaussian (gray circles in Fig. 2).

This proxy to the resolution is not ideal if the averaging kernel has a very complex shape, especially in case of smearing. For example, the gray circles on Fig. 2 Before choosing a value for η, we must also examine the influence on the model uncertainties.

Note that in this paper we express the uncertainties in the 1σ sense: there is a 68% probability or (i) would be appropriate if we were interested in smaller anomalies we expect to be strong enough to be detectable despite the increased uncertainties.

In the next subsection, we describe in more detail the η = 0.5 model, whose resolution lengths are on average short enough that the model is not too smooth, and whose uncertainties are on average low enough not to swamp all its features. Australia, where the resolution length is below 300 km; the worst resolved region away from the model's edges is located southeast of Hawaii, where the resolution length is greater than 1000 km.

Tomographic analysis with resolution and uncertainties 15

But the resolution is not only related to the path coverage: the data uncertainties also influence the achievable resolution and model uncertainties.

The uncertainty map in Fig. 3(e) shows that model uncertainties and model resolution do not always correlate. In some well-resolved regions, such as the Philippine Sea, Tonga-Kermadek, and central America, the locally-averaged shear-wave velocities are known to within 0.02 km/s. In these regions, we can resolve and interpret small features, even those that generate only moderate shear-wave velocity signatures. In some poorly-resolved regions, such as the central Pacific Ocean, shear-wave velocities are known to within greater bounds (0.07 km/s on average); here it becomes harder to interpret features unless they are large and generate strong shear-wave velocity signatures. Well-resolved regions may also display high model uncertainties: in the northern Pacific and off the coast of north-east Australia, the resolution length is close to 300 km, yet the uncertainties are similar to those in the central Pacific. In such regions, we can interpret features, regardless of their size, only if they generate strong shear-wave velocity signatures. This is a reminder that we need both resolution and uncertainty information for interpreting tomographic models.

The shear-wave velocity maps at 75 km depth (Fig. 3f), 125 km depth (Fig. 5c), and 400 km depth (Fig. 5f) show all the large-scale geological features we may expect in tomographic images in oceanic regions: low shear-wave velocities down to ∼100 km depth beneath the East Pacific Rise and the Pacific-Antarctic Ridge; fast velocities in the subduction zones (at least down to 200 km depth) and lower velocities in their back-arc regions; high velocities from 75 to ∼200 km depth in the North American and Australian cratons; velocities that increase with distance from the midocean ridges down to depths of ∼150 km. Unsurprisingly, given we have used their shear-wave velocity profiles, [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF] also saw similar features, though with slightly poorer resolution, as their damped-least-squares tomographic inversion used a single smoothing parameter for all locations in the model.

The resolution and uncertainty maps at 125 km and 400 km depth (Fig. 5) look like those at 75 km depth (Figs 3d ande). This is not what we would expect given the behaviour of surface waves. Deeper regions of the upper mantle can only be resolved by longer-period or higher-mode surface-waves, both of which have long wavelengths (so poor lateral resolution) and cannot be observed clearly on all paths; therefore, we would expect the resolution lengths to increase with depth. The reason those in Fig. 5 stay constant stems from the data we used in the SOLA inversions: path-averaged shear-wave velocity profiles from [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF], each one the result of a non-linear inversion of surface-wave measurements. Where long-period and higher-mode data were available, shear-wave velocity profiles were constrained by data down to 400 km depth;

where these data were unavailable, the profiles were damped, at depth, towards the background model: a smoothed version of PREM [START_REF] Dziewonski | Preliminary reference Earth model[END_REF]. It is likely, therefore, that many of the velocity profiles used for the SOLA inversion are uninformative below ∼300 km depth; we could not remove them from our inversion without redoing [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF]'s entire analysis, which was outside the scope of this study. In the following, therefore, bear in mind that the resolution maps (and for a similar reason the uncertainty maps) in Fig. 5 are overly optimistic at depths beyond ∼300 km.

Tomographic analysis with resolution and uncertainties 17

Comparing tomographic images with theoretical predictions

The tomographic maps of Fig 3 show an increase in seismic velocity with distance from the midoceanic ridges. This is connected with the well-known phenomenon of plate cooling (e.g. [START_REF] Ritzwoller | Cooling history of the Pacific lithosphere[END_REF][START_REF] Faul | The seismological signature of temperature and grain size variations in the upper mantle[END_REF][START_REF] Priestley | The thermal structure of the lithosphere from shear wave velocities[END_REF][START_REF] Isse | Surface wave tomography for the Pacific Ocean incorporating seafloor seismic observations and plate thermal evolution[END_REF]. Armed

with complete resolution and uncertainty information, we investigated if any features of these tomographic images deviated significantly from the predictions of a theoretical cooling model and if such deviations were well-resolved.

The reference model

We chose to use the simple half-space cooling model of [START_REF] Parker | Thermal Model of Ocean Ridges[END_REF] to illustrate our process of comparing tomographic images to theoretical predictions, not because we believed it to be the best thermodynamic fit to the ocean-cooling problem, but because it has the fewest adjustable parameters. Even with an unsuitable reference model, the analysis would still be relevant as it would simply lead us to state that the tomographic velocities differ significantly from the reference model.

The half-space cooling model describes the lithosphere as an infinite half-space of thermal diffusivity k that cools from a starting temperature Θ m . It predicts the temperature of the lithosphere as a function of age t and depth z from the following expression:

Θ(t, z) = Θ m erf z 1 4kt , ( 8
)
where we take Θ m = 1300 • C and k = 1.10 -6 m 2 s -1 , as suggested by [START_REF] Ritzwoller | Cooling history of the Pacific lithosphere[END_REF]. We turned this expression into temperature maps at each depth by relating age to geographic position using the age model of [START_REF] Müller | Age, spreading rates, and spreading asymmetry of the world's ocean crust[END_REF], then projected them onto our tomographic grid: the predicted temperature in cell j at depth z is given by Θ j (z) = Θ(t j , z), where t j is the average age of cell j.

We then had to relate temperatures predicted as a function of position with tomographic images that contained velocity estimates representing finite-resolution local averages. In order to make meaningful comparisons, we needed to account for this resolution and consider equivalent local averages of the predicted temperatures rather than the temperatures themselves. We obtained these local temperature averages using the resolution matrix :

Θ k (z) = M j=1 R (k) j (z)Θ j (z), a procedure
known as tomographic filtering (e.g. [START_REF] Ritsema | Tomographic filtering of geodynamic models: Implications for model interpretation and large-scale mantle structure[END_REF][START_REF] Simmons | Resolution and Covariance of the LLNL-G3D-JPS Global Seismic Tomography Model: Applications to Travel time Uncertainty and Tomographic Filtering of Geodynamic Models[END_REF].

In order to convert the locally-averaged predicted temperatures to shear-wave velocities, we needed to make some assumptions about the relationship between these two physical quantities that were at least partially supported by rock physics. We assumed that, in regions distant from Tomographic analysis with resolution and uncertainties 19 erence model exceeds the tomographic uncertainties by more than ±2σ m (panels (f) of the same figures). Since the estimation of the data uncertainties could be improved, so could the model uncertainties, we present ±1σ m and ±2σ m maps to appreciate a range of uncertainties.

Deviations from the reference model

However, just identifying anomalies as exceeding ±1σ m or ±2σ m is not enough to declare them significant, because even if the Earth were in reality identical to the prediction of the halfspace cooling model, we would still expect 32% of points in a tomographic model with uncertainties σ m to exceed ±1σ m and 5% of them to exceed ±2σ m . We could be justified in declaring anomalies to be significant only if more points than expected exceed the ±1σ m and ±2σ m thresh-Tomographic analysis with resolution and uncertainties 23 olds, or if these points organized geographically in coherent regions and these anomalous regions could indeed be resolved by the tomography (anomalies larger than the resolving lengths). This definition of significance is stricter than the one used in most tomographic studies, including [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF], and underlines the importance of correctly estimating the data uncertainties that feed into the estimates of σ m .

Table 2 shows the proportion of points that remain in the tomographic models after applying the ±1σ m and ±2σ m masks: at lithospheric depths, the proportions of such points do not exceed greatly the expected ones (0.32 and 0.05); however, as depth increases, so do the proportions of points exceeding the ±1σ m and ±2σ m levels. In section 2.1 we discussed that the data-uncertainties from [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF], and hence our estimates of σ m , may have been underestimated especially below ∼ 300 km. For this reason, we have limited our analysis of significant anomalies to depths shallower than ∼ 300 km.

Significant anomalies

At lithospheric depths (75 km), the proportion of points still visible in the tomographic image after applying the ±1σ m and ±2σ m masks is not significant (see Table 2). What about the geographic distribution of these points? Figs 6(e) and (f) shows that the positions of unmasked low-velocity anomalous regions at 75 km depth do not correlate with hot-spot locations published by [START_REF] Courtillot | Three distinct types of hotspots in the Earths mantle[END_REF], and that anomalies that exceed the ±2σ m thresholds are much smaller than the correlation lengths from Fig. 3(d). We therefore concluded that the lithosphere corresponding to the Pacific plate contains no significant shear-wave anomalies with respect to predictions made by the half-space cooling model, at least given the data set of [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF].

At shallow asthenospheric depths (200 km), a greater proportion of points remains in the masked tomographic image (see Table 2 and Figs 7e andf). A strong low velocity anomaly appears to the north-east of Hawaii at 20 • N 220 • E. This anomaly of 12 • (∼ 1400 km) radius is visible between 125 km and 250 km depth, and is about 0.25 km/s slower than the half-space cooling model (stronger than -3σ m ). Not only is the anomaly strong, it is also larger than the local resolution length of 750 km (Fig. 5a). Fig. 9 shows individual averaging kernels within and around the anomaly: the kernels are all of similar size or smaller than the anomaly and are all well focused, indicating absence of smearing. Given the amplitude and size of the anomaly with respect to the tomographic model's uncertainty σ m and its resolution length, we considered it to be significant.

We need to carefully analyse averaging kernels when interpreting the size and shape of anomalies.

Since the averaging kernels are spread in space, the velocity anomalies tend to appear larger than they really are. The shape of the anomaly in the tomographic image can also be misleading. For example the anomaly to the north-east of Hawaii shows two bumps, one to the north and one to the east, which may be due to the complex shape of the averaging kernels, or to data noise propagated into the model. Only features outside the uncertainty range with size greater than the resolution should be interpreted.

Deeper in the asthenosphere, a pattern of alternatively slow and fast velocity bands appears, oriented approximately south-east to north-west (Fig. 8). The bands are approximately 15 • wide (∼1700 km), 70 • (∼8000 km) long, and are visible between 275 km and 400 km depth. The bands seem to follow the absolute plate motion of the Pacific plate (e.g. [START_REF] Gripp | Current plate velocities relative to the hotspots incorporating the NUVEL-1 global plate motion model[END_REF]. They also resemble the bands observed by [START_REF] French | Waveform Tomography Reveals Channeled Flow at the Base of the Oceanic Asthenosphere[END_REF] at depths between 200 km and 350 km and the low anisotropy channels observed at 100 km depth by [START_REF] Montagner | Upper mantle low anisotropy channels below the Pacific Plate[END_REF]. The width and length of these bands are greater than the resolution length. To exclude possible smearing, always an issue when tomographic images show elongated features, we have shown several averaging kernels in Fig. 10. Given the amplitude and size of the highly-correlated bands and the focused nature of the averaging kernels, we again considered these bands to be significant features of the tomographic model.

DISCUSSION

We have produced a tomographic model of shear-wave velocities in the Pacific upper mantle using the SOLA Backus-Gilbert method [START_REF] Zaroli | Global seismic tomography using Backus-Gilbert inversion[END_REF][START_REF] Zaroli | Toward Seeing the Earth's Interior Through Unbiased Tomographic Lenses[END_REF], and have shown how to exploit the full model resolution and uncertainty information to evaluate the true significance of deviations from a theoretical prediction.

over the parameter space. Secondly these schemes drive the inversion to produce well-focused local averages (averaging kernels) that limit smearing along predominant paths (with SOLA this requires imposing well-focused target kernels, in our case circular ones). Thirdly, these schemes guarantee that the weights of the local averaging kernels (i.e. the rows of the resolution matrix) sum exactly to one, thereby eliminating the averaging bias that occurs with heterogeneous spatial distributions of data.

Model statistics -required for tomographic filtering and evaluating significance

Before comparing two different tomographic models, or a tomographic model to other data with different spatial resolution, we need to homogenize their resolutions, a process called tomographic filtering [START_REF] Ritsema | Tomographic filtering of geodynamic models: Implications for model interpretation and large-scale mantle structure[END_REF][START_REF] Simmons | Resolution and Covariance of the LLNL-G3D-JPS Global Seismic Tomography Model: Applications to Travel time Uncertainty and Tomographic Filtering of Geodynamic Models[END_REF]. This requires the knowledge of the full resolution matrix. In our study, we considered a simple temperature model of lithosphere cooling based on the half space cooling model of [START_REF] Parker | Thermal Model of Ocean Ridges[END_REF] that we filtered to the resolution of the tomographic model. Because the temperature model was very smooth, the tomographic filtering had only a minor effect. However, such filtering would be critical when comparing tomographic images with sharper models [START_REF] Mégnin | Imaging 3-D spherical convection models: What can seismic tomography tell us about mantle dynamics?[END_REF][START_REF] Simmons | Resolution and Covariance of the LLNL-G3D-JPS Global Seismic Tomography Model: Applications to Travel time Uncertainty and Tomographic Filtering of Geodynamic Models[END_REF][START_REF] Freissler | Tomographic filtering via the generalized inverse: a way to account for seismic data uncertainty[END_REF]).

To assess if some parts of a tomographic model significantly deviate from some reference we need to know the model uncertainties at each location. Because we had this information, we were able to mask deviations from the half-space cooling predictions that were smaller than ±1σ m and ±2σ m and argue for the significance of certain anomalies. The process we illustrated in section 3 could replace statistical appraisal techniques and other tomographic resolution tests that are not always well understood [START_REF] Lévêque | On the use of the checker-board test to assess the resolution of tomographic inversions[END_REF][START_REF] Rawlinson | Seismic Tomography and the Assessment of Uncertainty[END_REF][START_REF] Rawlinson | On the use of sensitivity tests in seismic tomography[END_REF].

Two significant anomalies in the Pacific region

We have focused on the significance of two anomalies: a low-velocity anomaly located to the north-east of Hawaii at 200 km depth and coherent bands of fast and slow velocities at 275 km depth.

Tomographic analysis with resolution and uncertainties 27

Anomaly NE of Hawaii

This emerges strongly from the background with an amplitude about 0.25 km/s slower than the average shear-wave velocity of the Pacific plate, where the uncertainty is around 0.07 km/s, thus corresponding to about -3σ m (Figs 7e and f). Despite using the same data-set, the model produced by [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF] shows a much weaker anomaly at the same location, probably because it was smoothed out by their damped-least-squares inversion and regularization. This has also been seen, though at weaker amplitudes than in this study, in tomographies constructed using fundamental mode Rayleigh waves with a finite-frequency scheme [START_REF] Liu | Global Rayleigh wave phase-velocity maps from finite-frequency tomography[END_REF] or multimode Love and Rayleigh waves [START_REF] Isse | Surface wave tomography for the Pacific Ocean incorporating seafloor seismic observations and plate thermal evolution[END_REF], while it was entirely missed by the GDM52 model of [START_REF] Ekström | A global model of Love and Rayleigh surface wave dispersion and anisotropy, 25-250 s: Global dispersion model GDM52[END_REF] based on fundamental mode Love and Rayleigh wave dispersion data.

Given the different data-sets, forward theories, inversion and damping schemes of these studies, it is hard to pinpoint the reason for their weaker anomaly, but we speculate that the anomaly was probably over-smoothed by the regularization of their inversions.

Coherent SE-NW bands

We see these alternating fast and slow bands emerging from the background at 275 km depth with deviations from HSC up to ±0.25 km/s (±1σ m and ±2σ m , see Fig. 8). [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF] do not show these bands in their 275 km depth image, suggesting that the SOLA inversion scheme was able to extract more information from the same dataset. Other studies have observed similar patterns: Ekström (2011) show a very weak pattern in their isotropic results, probably because fundamental-mode surface waves are only weakly sensitive at these depths; [START_REF] French | Waveform Tomography Reveals Channeled Flow at the Base of the Oceanic Asthenosphere[END_REF] see a stronger pattern with their full-waveform inversion; and, intriguingly, [START_REF] Liu | Global Rayleigh wave phase-velocity maps from finite-frequency tomography[END_REF] show similar bands in their dispersion maps at 30s and 50s, which are sensitive to much shallower depths.

Limits and perspectives

The SOLA inversion seems promising because it gives a way to control the model statistics (uncertainties and resolution) and to obtain them explicitly. The main problem is that it assumes

Figure 1 .

 1 Figure 1. Data geometry and target resolution at 75 km depth. Thick green lines indicate plate boundaries. (a) Each thin line represents one of the ∼15 000 paths from Maggi et al. (2006) used as input to the SOLA tomography; the color of each path is the path-average shear wave-speed at 75 km depth in the 1D model corresponding to the path; (b) Path density (total path length that falls in each cell normalized by the total path length over the entire model); (c) Target-kernel radii.

  Figs 1(b) and (c) show the path densities and target resolution length calculated for the path distribution in Fig. 1(a). We show several target kernels in Fig. 2(a).

3. 1

 1 Resolution, uncertainty, and model estimate for various trade-off parametersAll tomographic inversions have trade-offs. In damped-least-squares inversions, model smoothness trades-off against data fit (smooth models fit the data poorly, but rough models risk fitting the noise part in the data). In SOLA, model resolution trades-off against model uncertainty. Fig.2

Figure 2 .

 2 Figure 2. Examples of SOLA target kernels T (k) (a) and corresponding averaging kernels A (k) (b-d) at 75 km depth for η=3, 0.5 and 0.01. Blue circles indicate the edge of the target kernel disks; grey circles contain 68% of the averaging kernel amplitude. For convenience, several target (or averaging) kernels are plotted on the same map.

  (b) do not represent well the averaging kernels when the ramifications are long and strong (especially for those in the middle of the Pacific plate). Note that other simplifications of the resolution could be used; for instance we could fit an ellipse to the averaging kernels to obtain a main direction in case of anisotropic resolution (smearing). However, since our averaging kernels are mostly circular, we use the circular approximation for simplicity. The top row of Fig.3shows resolution maps for three values of the trade-off parameter. In all maps the resolution is bad at the model borders and in the central Pacific Ocean, where path coverage is poor. As the trade-off parameter decreases, resolution improves in the central Pacific, though the improvement seems to level off between η = 0.5 and η = 0.01.

Figure 3 .

 3 Figure 3. SOLA inversion results at 75 km depth for three values of trade-off parameters: η=3, 0.5, 0.01. Each column corresponds to a value of η and from top to bottom, maps show resolution lengths (L k ) (a, d, g), model uncertainties (b, e, h) and shear-wave velocities (c, f, i).

Figure 4 .

 4 Figure 4. Mean resolution versus mean uncertainty of SOLA results for five values of the trade-off parameter η (5, 3, 1, 0.5, 0.01) indicated on the curve. The slope from η = 3 to η = 0.5 is more negative than from η = 0.5 to η = 0.01. This is consistent with the apparent variations of the model statistics on Fig. 3. As η changes, model resolution and uncertainties follow an L-curve.

Fig. 4

 4 Fig. 4 shows mean resolution ( 1 M

3. 2

 2 Detailed analysis of the resolution, uncertainty and model estimate at a given trade-off For any tomographic inversion, resolution, model uncertainty, and model smoothness are intrinsically linked to path coverage and data quality. The best-sampled regions for our inversion are the eastern, northern and western borders of the Pacific ocean (Fig. 1): these regions are strongly seismogenic and host many seismic stations. The interior of the Pacific ocean and its southern boundary are sparsely sampled, except near Hawaii and French Polynesia where Maggi et al. (2006) exploited both permanent stations and temporary deployments (Barruol 2002) to increase local path coverage. The resolution map in Fig. 3(d) shows that the well-sampled regions have the best lateral resolution: the radii of the circles that approximate the averaging kernels in these regions are between 300 km and 500 km. The Pacific ocean itself has a poor resolution, around 800 km on average. The regions west of Hawaii and west of the French Polynesia have better resolution (550 km) because the seismic networks installed within those archipelagos record earthquakes that occur on the western boundary of the Pacific plate. The best resolved region is northeast

Figure 5 .

 5 Figure 5. SOLA inversion results at 125 km (first row) and 400 km depth (second-row) for η = 0.5: (a, d) resolution lengths; (b, e) model uncertainties; (c, f) shear-wave velocities.

Figure 6 .

 6 Figure 6. Comparison at 75 km depth between the SOLA tomographic model and the predictions of the half-space cooling (HSC) model from Parker & Oldenburg (1973) on the Pacific plate, excluding ridges and subduction zones. a) shear-wave velocities from SOLA; b) shear-wave velocities predicted for the HSC model; c) difference between SOLA and HSC; d) difference between SOLA and HSC scaled by the tomographic uncertainties σ m ; e) and f) same as d) but with masks to remove deviations smaller than ±1σ m or ±2σ m respectively. Grey contour lines indicate lithosphere age from Müller et al. (2008) and yellow triangles locate hot-spots from Courtillot et al. (2003).

Figure 7 .

 7 Figure 7. Same as Fig. 6 but at 200 km depth (note that at this depth the HSC model shows no agedependence). The black line indicates the contour of the anomaly located to the north-east of Hawaii.

Figure 8 .

 8 Figure 8. Same as Fig. 6 but at 275 km depth. The black lines indicate low velocity anomalous bands (the lines follow approximately the zeroes between anomalous bands).

Figs 6

 6 Figs 6(c), 7(c), and 8(c) show the differences between the tomographic model and the HSC reference model at 75 km, 200 km, and 275 km depth. Red anomalies correspond to areas where the tomographic model is slower than HSC; blue anomalies correspond to areas where the to-

Figure 9 .

 9 Figure 9. Selection of averaging kernels around the slow velocity anomaly located to the north-east of Hawaii (see Fig. 7) at 200 km depth. The black line represents the anomaly contour.

Figure 10 .

 10 Figure 10. Selection of averaging kernels in the Pacific at 275 km depth. The black lines represent the pattern of anomalous velocity bands observed on Fig. 8. The averaging kernels are located at: ( 30N, 165E), ( 30N, 195E), ( 30N, 225E), ( 0N, 165E), ( 0N, 195E), ( 0N, 225E), ( -30N, 165E), ( -30N, 195E), ( -30N,225E). For convenience, we plot the averaging kernels on the same map.

Table 1 .

 1 χ 2 red values at various depths for the damped-least-squares model of[START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF].

	depth [km]	50	75	100	150	200	250	300	350	400
	χ 2 red	3.378 3.199 2.630 2.216 1.648 1.561 1.542 1.579 1.497
	mates the full data uncertainties. Because the data uncertainties influence the SOLA inversion and
	propagate directly into model uncertainties (see section 2.2), we needed more reliable estimates of
	them. We calculated the χ 2 red for the final model of Maggi et al. (2006) (e.g. Nolet 2008):

Table 2 .

 2 Proportion of points remaining after masking values within the 1σ m or 2σ m error bar(s) at various depths. For a normal distribution, we would expect proportions of 0.32 for 1σ m and 0.05 for 2σ m .

	Depth (km) 50	75	100 150 200 250 300 350 400
	For 1σ

m 0.32 0.33 0.33 0.33 0.38 0.41 0.46 0.49 0.43 For 2σ m 0.06 0.03 0.03 0.07 0.09 0.10 0.14 0.13 0.12

Acknowledgements

The authors are grateful to Huajian Yao (editor), Frederik Tilmann, and an anonymous reviewer whose constructive comments helped to improve this paper.

DATA AND CODE AVAILABILITY

The computer codes for the SOLA inversions and the data used in this study are available from C.Z.

(c.zaroli@unistra.fr) and A.M. (alessia.maggi@unistra.fr), respectively, upon reasonable request.

Several authors have inverted surface-wave data to obtain velocity models of the upper mantle, either globally or in the Pacific, using different forward theories and inversion methods (e.g. [START_REF] Zhou | Global upper-mantle structure from finite-frequency surface-wave tomography[END_REF][START_REF] Ekström | A global model of Love and Rayleigh surface wave dispersion and anisotropy, 25-250 s: Global dispersion model GDM52[END_REF][START_REF] French | Waveform Tomography Reveals Channeled Flow at the Base of the Oceanic Asthenosphere[END_REF][START_REF] Liu | Global Rayleigh wave phase-velocity maps from finite-frequency tomography[END_REF][START_REF] Isse | Surface wave tomography for the Pacific Ocean incorporating seafloor seismic observations and plate thermal evolution[END_REF]. Although they agree with each other at large scales, they often display different local features. These discrepancies may be partly explained by the differing data-sets (fundamental-mode Rayleigh waves, multi-mode Rayleigh waves, or full-waveforms) and inversion schemes, but may also be due to the complex relationship between the heterogeneous path coverage of the Pacific region, irregular resolution, and data-errors.

Without full resolution and uncertainty information -model statistics -we cannot give strong quantitative arguments for or against the significance of any small-scale feature of a tomographic model (e.g. [START_REF] Foulger | Caveats on tomographic images[END_REF]. Unfortunately, most damped-least-squares or iterative conjugategradient inversion schemes do not provide model statistics or control them during the inversion, as they focus instead on minimising the misfit between observed measurements and those predicted by the tomographic model. It is possible to approximate the model statistics for these inversion schemes (see for example [START_REF] Rawlinson | Seismic Tomography and the Assessment of Uncertainty[END_REF][START_REF] Rawlinson | On the use of sensitivity tests in seismic tomography[END_REF], for reviews of various techniques), but such methods are rarely used in practice as they are computationally costly and often provide only crude estimates of the model statistics. For the tomography problem considered here, it would indeed be possible to explicitly compute the G † matrix from a 'traditional' damped-least-squares inversion and then propagate the data uncertainties into model uncertainties; but the advantage of the SOLA method is that it allows a direct control on the model resolution and uncertainties.

Advantages of Backus-Gilbert type inversions

Beyond delivering model resolution and uncertainties at no extra cost, there are other important advantages of using inversion schemes based directly on the model statistics, such as those proposed by [START_REF] Backus | Numerical Applications of a Formalism for Geophysical Inverse Problems[END_REF][START_REF] Backus | The Resolving Power of Gross Earth Data[END_REF][START_REF] Backus | Uniqueness in the inversion of inaccurate gross Earth data[END_REF], [START_REF] Pijpers | Faster formulations of the optimally localized averages method for helioseismic inversions[END_REF], 1993), [START_REF] Zaroli | Global seismic tomography using Backus-Gilbert inversion[END_REF], and [START_REF] Zaroli | Toward Seeing the Earth's Interior Through Unbiased Tomographic Lenses[END_REF]. Firstly, these inversion schemes constrain model uncertainties to vary smoothly Latallerie et al.

the forward problem is linear while surface-wave physics is not (the dispersion characteristics of surface-waves are non-linearly related to the shear-wave velocity as a function of depth). Surfacewave tomography is often performed in two steps: a non-linear step to relate frequency-dependent phase and/or group velocities to depth-dependent shear-wave velocity; and a linear regionalization step. Depending on the order of the two steps, the regionalization is performed either on group and phase velocities, or on 1D shear-wave velocities as in [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF]. Only the linear step can be performed using SOLA.

In this study, as in [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF], each depth in the tomographic model was inverted independently, so vertical coherence could not be imposed or evaluated using resolution kernels.

Furthermore, as the inputs to the SOLA inversion were 1D path-averaged shear-wave velocity models defined on the same depth range, the geographical coverage was identical for each depth, and did not take into account the greater sensitivity to depth of some paths (illuminated by longerperiod or higher-mode surface waves) with respect to others.

To obtain a fully 3D shear-wave velocity model directly from frequency-dependent measurements of surface-waves in a single linear step with SOLA, and to take the full depth-sensitivity of individual measurements into account, we would need to adopt a finite-frequency description of surface waves to help linearize the tomographic problem (e.g. [START_REF] Zhou | Three-dimensional sensitivity kernels for surface wave observables[END_REF][START_REF] Yoshizawa | Sensitivity kernels for finite-frequency surface waves[END_REF]. This would also allow us to obtain a fully 3D model with 3D resolution kernels, and therefore be able to apply the same statistical rigor to interpreting any vertical structures within the model (e.g. plumes, the lithosphere depth).

Like any other process, the pertinence of the SOLA outputs rest upon the quality of the inputs.

In particular, the model errors we have relied upon to decide the significance of anomalies are nothing more than a propagation into the model space of the data errors. While traditional dampedleast-squares schemes use data uncertainties to weigh data relatively to each other and to evaluate model quality (using χ 2 red measures, for example), SOLA uses them more directly to evaluate uncertainties in the model which then drive the inversion through their trade-off with resolution.

In order to be able to trust the magnitude of the model errors in SOLA, it is necessary for the magnitude of the data errors to be correct. This is the reason why we upscaled the data uncertainties Tomographic analysis with resolution and uncertainties 29 using the χ 2 red values of the damped-least-squares model of [START_REF] Maggi | Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature[END_REF]; however a more rigorous assessment of the data uncertainties is still required.

Our study calls for other surface wave tomography studies based on the SOLA inversion in other regions of the world, with new data sets, and particular care in estimating data uncertainties.

Our tomographic model, together with its resolution and uncertainty maps, calls for new data in the Pacific to observe other seismic velocity anomalies significant enough to be interpreted robustly.

CONCLUSION

Rayleigh wave tomography of the Pacific upper mantle shows large scale geological features (e.g. lithosphere cooling) but both the model uncertainties and resolution are required for interpretations at smaller scales. We exploited the SOLA tomographic inversion scheme [START_REF] Zaroli | Global seismic tomography using Backus-Gilbert inversion[END_REF] to propose a workflow to analyse the tomographic model using its resolution and uncertainties, that is:

(i) Define some reference, Earth-like, physical model;

(ii) Filter the reference model to the tomographic resolution;

(iii) Compute the deviations from the reference model;

(iv) Normalize the deviations by the model uncertainties;

(v) Mask model estimates that are within one or two error bars to focus only on significant seismic anomalies;

(vi) Compare those non-masked anomalies to the local resolution to discuss their resolvability and spot artifacts, if any.

In this study, we used a half space cooling description of the Pacific lithosphere as a reference model. We observed a low velocity anomaly to the north-east of Hawaii at 200 km depth and a pattern of alternatively slow and fast anomalous bands oriented approximately south-east to northwest at 275 km depth. Both features are reliable given the model statistics. Our study shows that the model statistics can be properly analysed in surface wave tomography using the SOLA inversion;

the same scheme could be applied in other regions of the globe.

In this study, the data uncertainties were estimated using multiple sampling with ray paths close to each other, then upscaled because apparently underestimated. We should aim collectively to characterize more accurately uncertainties in seismological data sets -clearly a challenging but crucial task.


Niklas Wellander

Sebastien Guenneau

Elena Cherkaev

TWO-SCALE CUT-AND-PROJECTION CONVERGENCE FOR QUASIPERIODIC MONOTONE OPERATORS

Averaging certain class of quasiperiodic monotone operators can be simplified to the periodic homogenization setting by mapping the original quasiperiodic structure onto a periodic structure in a higher dimensional space using cut-and projection method. We characterize cut-and-projection convergence limit of the nonlinear monotone partial differential operator -div σ x, Rx η , ∇uη for a bounded sequence uη in W 1,p 0 (Ω), where 1 < p < ∞, Ω is a bounded open subset in R n with Lipschitz boundary. We identify the homogenized problem with a local equation defined on the hyperplane in the higher-dimensional space. A new corrector result is established.

Introduction

Nonlinear physical phenomena are ubiquitous in modern electronic devices. A few examples are current surge protectors made of varistor ceramics, solid state amplifiers and integrated circuits. This is one motivation to develop mathematical tools that can be used to analyse the effective properties of polycrystalline quasiperiodic semiconductors.

In [START_REF] Braides | Homogenization of penrose tilings[END_REF], it is shown that integral energies F η where the spatial dependence follows the geometry of a Penrose tiling, or more general quasicrystalline geometries, can be homogenized. More precisely,

(1)

F η (x) = Ω f x η , ∇u(x) dx , u ∈ W 1,p (Ω)
where Ω is an open subset of R 2 , and f depends on x through the shape and the orientation of the cell containing x in an a-periodic tiling of the space, Γ-converge on W 1,p (Ω) with respect to the L p convergence to the functional

(2) F 0 (x) = lim inf

T →∞ 1 T 2 (0,T ) 2 f (y, ∇v(y) + ξ) dy , v ∈ W 1,p 0 ((0, T ) 2 )
This general homogenization result was shown using that f is Besicovitch almost periodic in y and thus a previous result on Besicovitch almost periodic functionals [START_REF] Braides | A homogenization theorem for weakly almost periodic functionals[END_REF] could be applied. Homogenization of interfacial energies on Penrose lattices making use of Γ-convergence for similar functionals to (1) but with the volume integral replaced by a surface integral has been also addressed in [START_REF] Braides | Interfacial energies on quasicrystals[END_REF]. Γ-convergence is a very powerful tool in homogenization theory [B + 02], but twoscale convergence [START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF][START_REF] Allaire | Homogenization and two-scale convergence[END_REF] can more easily identify homogenized equations in the periodic setting. A similar tool is the periodic unfolding approach [START_REF] Cioranescu | The periodic unfolding method in homogenization[END_REF], in which one first maps the original sequence of functions to a sequence that is defined on R n ×]0, 1[ n , and then takes the usual weak limit in suitable function spaces, using this extended domain. This is similar to the two-scale Fourier transform approach proposed in [START_REF] Wellander | The two-scale Fourier transform approach to homogenization; periodic homogenization in Fourier space[END_REF].

Due to its simplicity, one might wish to apply two-scale homogenization (or the periodic unfolding and Fourier transform approaches) to quasiperiodic materials or mixtures of materials with rational and irrational periodicity, e.g., see [START_REF] Braides | Correctors for the homogenization of almost periodic monotone operators[END_REF][START_REF] Casado-Diaz | The two-scale convergence method applied to generalized besicovitch spaces[END_REF] for a setting in an almost periodic regime [START_REF] Blanc | Local profiles for elliptic problems at different scales: defects in, and interfaces between periodic structures[END_REF][START_REF] Blanc | Profils locaux et problèmes elliptiques à plusieurs échelles avec défauts[END_REF] for some recent work on quasiperiodic multiscale homogenization setting.

This has been proposed in [START_REF] Guenneau | Homogenisation of quasi-crystals and analysis of modes in photonic crystal fibres[END_REF][START_REF] Bouchitté | Homogenization of dielectric photonic quasi crystals[END_REF] wherein two-scale convergence is applied to the quasi-periodic setting making use of the cut-and-projection method. Indeed, quasiperiodic materials can be often described by periodic structures in higher spatial dimensions cut and projected onto a hyper-plane or a lower dimensional (physical) space, typically R 2 (such as for Penrose tilings) and R 3 , as proposed by the physicists Duneau and Katz [START_REF] Duneau | Quasiperiodic patterns[END_REF]. This makes it possible to use standard periodic homogenization tools such as two-scale convergence, to homogenize quasiperiodic materials.

To do that one has to complement existing tools with the cut-and-projection operator, this was done in [START_REF] Bouchitté | Homogenization of dielectric photonic quasi crystals[END_REF] in the framework of W 1,2 , making use of Fourier representation for two-scale limits of gradients. Importantly, this was revisited in [START_REF] Wellander | Two-scale cut-and-projection convergence; homogenization of quasiperiodic structures[END_REF][START_REF] Wellander | Homogenization of quasiperiodic structures and two-scale cut-and-projection convergence[END_REF].

However, two-scale convergence can also be applied to nonlinear operators [START_REF] Allaire | Homogenization and two-scale convergence[END_REF]. This is a generalization of the usual weak convergence in Lebesgue spaces L p , 1 < p < ∞, in which one uses oscillating test functions to capture oscillations on the same scale as the test functions in the sequence of functions that are investigated. As a consequence one obtains limit functions that are defined on the product space

R n ×]0, 1[ n .
In this paper we extend the two-scale cut-and-projection convergence method to Sobolev spaces W 1,p , 1 < p < ∞. We build upon [START_REF] Wellander | Homogenization of quasiperiodic structures and two-scale cut-and-projection convergence[END_REF] to characterize the limits of nonlinear partial differential operators in this setting. We illustrate the method on a nonlinear electrostatic problem that was previously homogenized using the tool of G-convergence for a larger class of almost periodic functions in [START_REF] Braides | Homogenization of almost periodic monotone operators[END_REF]. We finally establish a corrector result for the gradients. We note that Kozlov established in [START_REF] Mikhailovich | Averaging differential operators with almost periodic, rapidly oscillating coefficients[END_REF] a corrector result that requires almost periodic coefficients, which are the restriction of sufficiently smooth periodic functions of greater number of variables, when the problem is set in the whole space.

1.1. Setup of a nonlinear electrostatic equation in a quasiperiodic structure. Throughout the paper, we consider a bounded domain Ω in R n with Lipschitz boundary. We study the electrostatic equation, which is applicable to model DC currents in semiconductors, e.g., ZnO based varistor ceramics,

(3)

   -div σ x, Rx η , ∇u η (x) = f (x) , x ∈ Ω u η | ∂Ω = 0
with f ∈ W -1,q (Ω) and p is a real constant 1 < p < ∞, and q is its dual exponent, i.e., 1/p + 1/q = 1. We use standard notations for Lebesgue and Sobolev spaces. The Euclidean norm and the scalar product in R n are denoted by | • | and (•, •), respectively. The current density is given by a non-linear map, σ, that satisfies assumptions i)-vii

): i) σ(x, •, ξ) is Y -periodic in R m and Lebesgue measurable for every x, ξ ∈ R n . ii) σ(•, y, ξ) is continuous for almost every y ∈ R m and every ξ ∈ R n . iii) σ(x, y, •) is continuous for almost every x ∈ R n and y ∈ R m . iv) 0 ≤ c|ξ| p ≤ (σ(x, y, ξ), ξ) , c > 0 for almost every x ∈ R n and y ∈ R m . v) (σ(x, y, ξ 1 ) -σ(x, y, ξ 2 ), ξ 1 -ξ 2 ) ≥ c 1 |ξ 1 -ξ 2 | p , c 1 > 0 for all ξ 1 , ξ 2 ∈ R n , and almost every x ∈ Ω and y ∈ R m . vi) |σ(x, y, ξ)| ≤ c 2 1 + |ξ| p-1 , c 2 > 0, ∀ξ ∈ R n , almost every x ∈ Ω and y ∈ R m . vii) σ(•, R • /η, ξ) is Lebesgue measurable, for all ξ ∈ R n .
Standard estimates yield solutions that are uniformly bounded in W 1,p 0 (Ω) with respect to η.

Remark 1.1. Assumption v) is needed for a corrector result. It can be replaced by a strict monotone assumption if only the homogenized equation is needed, e.g., see [START_REF] Allaire | Homogenization and two-scale convergence[END_REF].

1.2. Two-scale cut-and-projection convergence. In this section, we recall some properties of two-scale convergence in L p (Ω), 1 < p < ∞, Ω ⊂ R n [START_REF] Allaire | Homogenization and two-scale convergence[END_REF], and revisit the extension to the quasiperiodic setting, see [START_REF] Bouchitté | Homogenization of dielectric photonic quasi crystals[END_REF] when p = 2. More precisely, we consider a real valued matrix R with m rows and n columns, and we would like to approximate an oscillating sequence {u η (x)} η∈]0,1[ by a sequence of two-scale functions u 0 x, Rx η where u 0 (x, •) is Y m -periodic on R m , a higher dimensional space paved with periodic cells

Y m = ]0, 1[ m .
In what follows, we assume that R : R n → R m , n < m, fulfills the criterion

(4) R T k = 0 , ∀k ∈ Z m \ {0}
In fact, when defining a quasiperiodic structure through cut-and-projection, one notes that the matrix R is not uniquely defined (e.g., an icosahedral phase using a mapping from R 6 or R 12 onto R 3 , or a Penrose tiling using a mapping from R 4 or R 5 onto R 2 , see [START_REF] Duneau | Quasiperiodic patterns[END_REF][START_REF] Janot | Quasicrystals: A Primer[END_REF]). However, if g is a trigonometric polynomial, then the f = g • R admits the following (uniquely defined) ergodic mean:

(5)

L(f ) = lim T →+∞ 1 (2T ) n ]-T ;T [ n f (x) dx = Y m g(y) dy = [g]
where [g] denotes the mean of g over the periodic cell Y m in R m . This is the case provided that R fulfills the criterion (4), see [START_REF] Bouchitté | Homogenization of dielectric photonic quasi crystals[END_REF]. This result suggests the following concept of two-scale convergence attached to a matrix R. Definition 1.1 (Distributional two-scale convergence). We say that the sequence

(u η ) in L p (Ω), 1 < p < ∞, two-scale converges in the distributional sense towards the function u 0 ∈ L p (Ω × Y m ) for a matrix R, if for every ϕ ∈ D(Ω; C ∞ (Y m )): (6) lim η→0 Ω u η (x)ϕ x, Rx η dx =
Ω×Y m u 0 (x, y)ϕ(x, y) dxdy Definition 1.2 (Weak two-scale convergence). We say that the sequence (u η ) in L p (Ω), two-scale converges weakly towards the function

u 0 ∈ L p (Ω × Y m ) for a matrix R, if for every ϕ ∈ L q (Ω, C (Y m )), 1 < p < ∞, 1/p + 1/q = 1, (6) holds.
We denote weak two-scale convergence for a matrix R with u η k R u 0 . The following result, which is a straightforward extension of a proof in [START_REF] Bouchitté | Homogenization of dielectric photonic quasi crystals[END_REF] to L p case corresponding to Corollary 1.15 in [START_REF] Allaire | Homogenization and two-scale convergence[END_REF], ensures the existence of such twoscale limits when the sequence (u η ) is bounded in L p (Ω) and R satisfies (4).

Proposition 1.1. If R is a matrix satisfying ( 4) and (u η ) is a bounded sequence in L p (Ω), 1 < p < ∞, then there exist a vanishing subsequence η k and a limit

u 0 (x, y) ∈ L p (Ω × Y m ) (Y m -periodic in y) such that u η k R u 0 , as η k → 0.
A proof of Proposition 1.1 uses the same arguments as in the periodic case, e.g., see [START_REF] Lukkassen | Two-scale convergence[END_REF] and can be found in [START_REF] Ferreira | Homogenization of quasi-crystalline functionals via two-scale-cut-and-project convergence[END_REF].

We will need to pass to the limit in integrals

Ω u η v η dx where u η R u 0 and v η R v 0 .
For this, we introduce the notion of strong two-scale (cut-and-projection) convergence for a matrix R.

Definition 1.3 (Strong two-scale convergence).

A sequence u η in L p (Ω) is said to two-scale converge strongly, for a matrix R, towards a limit u 0 in L p (Ω × Y m ),

which we denote u η R → → u 0 , if and only if u η k R u 0 and (7) u η (x) L p (Ω) → u 0 (x, y) L p (Ω×Y m )
This definition expresses that the effective oscillations of the sequence (u η ) are on the order of η. Moreover, these oscillations are fully identified by u 0 .

We note the following result which is useful to establish a link between weak quasiperiodic convergence in definition 1.2 and strong L p convergence (see Corollary 1.2).

Proposition 1.2. Let R be a linear map from R n to R m satisfying (4). Let f (x, y) ∈ L 1 (Ω; C (Y m )). Then f x, Rx η is a measurable function on Ω such that (8) f x, Rx η L 1 (Ω) ≤ f (x, y) L 1 (Ω;C (Y m )) := Ω sup y∈Y m | f (x, y) | dx and (9) lim η→0 Ω f x, Rx η dx = Ω×Y m f (x, y) dxdy Proof Measurability of f x, Rx η follows from Theorem 1 in [LNW02] that
ensures f is of Carathedory type. Inequality (8) is obvious. Finally, (9) follows from Lemma 2.4 in [START_REF] Bouchitté | Homogenization of dielectric photonic quasi crystals[END_REF].

We deduce two corollaries. First, we have a result that is useful to establish a corrector type result (see proposition 1.3).

Corollary 1.1. Let R be a linear map from R n to R m satisfying (4). Let φ(x, y) ∈ L p (Ω; C (Y m )). Then (10) lim η→0 Ω φ x, Rx η p dx = Ω×Y m |φ(x, y)| p dxdy
Proof We first consider the case when φ can be expressed as

φ(x, y) = τ (x)β(y), where τ (x) ∈ L ∞ (Ω) and β(y) ∈ C (Y m ). Since β p ∈ C (Y m ), we deduce from proposition 1.2 that β p Rx η converges towards its mean [β p ] weakly in L 1 (Ω). As τ p belongs to L ∞ (Ω), we obtain (11) lim η→0 Ω τ p (x)β p Rx η dx = Ω τ p (x) dx Y m β p (y)dy
From Fubini's theorem, this implies that (10) holds. This result is extended by linearity to step functions

φ k ∈ S t (Ω, C (Y m )) such that φ k = k i=1 t i χ Ai (x)ψ i (y), where A i = {x ∈ Ω , φ k (x, .) = t i } and ψ i (y) ∈ C (Y m
). We deduce that (10) holds by density in L p (Ω, C (Y m )).

More precisely, we consider φ k ∈ S t (Ω, C (Y m )). There exists a sequence of step

functions φ k = k i=1 t i χ Ai (x)ψ i (y) such that (12) lim k→∞ Ω sup y∈Y m | φ k (x, y) -φ(x, y) | p dx = lim k→∞ φ k -φ p L p (Ω,C (Y m ))
Moreover, from the triangular inequality and the continuity of the linear map R, we deduce that there exists a constant C > 0 such that (13)

φ x, Rx η L p (Ω) ≤ C φ x, Rx η -φ k x, Rx η L p (Ω) + φ k x, Rx η L p (Ω)
Noting that for every

v ∈ L p (Ω, C (Y m )) (14) Ω v x, Rx η p dx ≤ Ω sup y∈Y m |v(x, y)| p dx = v p L p (Ω,C (Y m ))
we deduce from (13) that for every integer k

(15) φ x, Rx η L p (Ω) ≤ C φ -φ k L p (Ω,C (Y m )) + φ k x, Rx η L p (Ω)
Since φ k is admissible, we deduce that there exists a constant C > 0 such that (16)

lim sup η→0 φ x, Rx η L p (Ω) ≤ C φ -φ k L p (Ω,C (Y m )) + φ k L p (Ω×Y m ) ≤ C φ -φ k L p (Ω,C (Y m )) + φ k -φ L p (Ω×Y m ) + φ L p (Ω×Y m ) ≤ (C + 1) φ -φ k L p (Ω,C (Y m )) + φ L p (Ω×Y m )
Passing to the limit k → ∞ in (16), we obtain

(17) lim sup η→0 φ x, Rx η L p (Ω) ≤ φ L p (Ω×Y m )
Similar arguments hold for the lower limit. In that case, we obtain that for every integer k (18)

lim inf η→0 φ x, Rx η L p (Ω) ≥ φ k L p (Ω×Y m ) -C φ -φ k L p (Ω,C (Y m )) ≥ φ L p (Ω×Y m ) -φ k -φ L p (Ω×Y m ) -C φ -φ k L p (Ω,C (Y m )) ≥ φ L p (Ω×Y m ) -(C + 1) φ -φ k L p (Ω,C (Y m ))
Passing to the limit k → ∞ in (18), we obtain

(19) lim inf η→0 φ x, Rx η L p (Ω) ≥ φ L p (Ω×Y m )
Equality ( 10) is established combining inequalities ( 17) and (19).

A second Corollary of proposition 1.2 is Corollary 1.2. Let Ω be an open bounded set in R n and Y m = ]0, 1[ m with m > n. Let a sequence (u η ) converge strongly to u in L p (Ω). Then there exist a vanishing subsequence η k and a limit u 0 (x, y)

∈ L p (Ω × Y m ) (Y m -periodic in y) such that u η k R u 0 = u as η k → 0 for any m × n matrix R satisfying (4). Proof Let φ(x, y) ∈ L q (Ω; C (Y m )). Then, we have (20) Ω u η (x)φ x, Rx η dx - Ω×Y m u 0 (x, y)φ(x, y) dxdy ≤ u η (x) -u(x) L p (Ω) φ x, Rx η L q (Ω) + Ω u(x)φ x, Rx η dx - Ω×Y m u(x)φ(x, y) dxdy Since u(x)φ(x, y) ∈ L 1 (Ω; C (Y m )) as u is in L p (Ω)
by assumption, proposition 1.2 ensures the second term in RHS goes to zero. The first term in RHS goes to zero since by corollary 1.1

(21) lim η→0 Ω φ x, Rx η q dx = Ω×Y m | φ(x, y) | q dxdy
and (u η ) converges strongly to u in L p (Ω) by assumption.

The following proposition provides us with a corrector type result for the sequence u η when its limit u 0 is smooth enough:

Proposition 1.3. Let R be a linear map from R n in R m satisfying (4). Let u η be a sequence such that u η R u 0 (x, y) (weakly). Then i) u η weakly converges in L p (Ω) towards u(x) = Y m
u 0 (x, y) dy and

(22) lim inf η→0 u η L p (Ω) ≥ u 0 L p (Ω×Y m ) ≥ u L p (Ω)
ii) Let v η be another bounded sequence in L q (Ω), 1/p + 1/q = 1, such that v η R → → v 0 (strongly), then

(23) u η v η → w(x) in D (Ω) where w(x) = Y m u 0 (x, y)v 0 (x, y) dy iii) If u η R → → u 0 (x, y) strongly and u 0 is smooth enough (e.g., u 0 ∈ L p (Ω, C (Y m ))) and (24) u 0 x, Rx η L p (Ω) → u 0 (x, y) L p (Ω×Y m ) then (25) u η -u 0 x, Rx η L p (Ω) → 0
Proof (i) Choosing test functions φ in L q (Ω; C (Y m )) independent of the y variable in definition 1.2, one has that for every φ ∈ L q (Ω)

(26) lim η→0 Ω u η (x)φ (x) dx = Ω×Y m u 0 (x, y)φ(x) dxdy = Ω φ (x)
Y m u 0 (x, y)dy dx Moreover, (u η ) is bounded in L p (Ω) as a weakly convergent sequence in a Hilbert space.

Then, let ϕ m be a sequence in

L q (Ω, C (Y m )) such that ϕ m converges to | u 0 | p-2 u 0 strongly in L q (Ω × Y m ).
We first apply the Young inequality for real numbers a and b, and 1 < p < ∞,

1/p + 1/q = 1, which states that ab ≤| a | p /p+ | b | q /q. We consider a = u η and b = ϕ m to get (27) Ω | u η (x) | p dx ≥ p Ω u η (x)ϕ m x, Rx η dx -(p -1) Ω ϕ m x, Rx η q dx
We first pass to the limit when η goes to zero: (28)

lim inf η→0 u η p L p (Ω) ≥ p Ω×Y m u 0 (x, y)ϕ m (x, y) dxdy-(p-1) Ω×Y m |ϕ m (x, y)| q dxdy
where we have used that u η R u 0 (x, y) weakly and

ϕ m x, Rx η R → → ϕ m (x, y)
strongly (making use of proposition 1.1). We then pass to the limit when m goes to infinity:

(29) lim inf η→0 u η p L p (Ω) ≥ p Ω×Y m | u 0 (x, y) | p dxdy-(p-1) Ω×Y m | u 0 (x, y) | p dxdy = u 0 p L p (Ω×Y m )
where we have used that

ϕ m →| u 0 | p-2 u 0 strongly in L q (Ω × Y m ).
Moreover, thanks to Jensen's inequality, we have that (30)

u p L p (Ω) = Ω Y m u 0 (x, y)dy p dx ≤ Ω×Y m | u 0 (x, y) | p dxdy = u 0 p L p (Ω×Y m )
We conclude by combining (29) and ( 30).

(ii) Let ψ m be a sequence in

L p (Ω, C (Y m )) such that ψ m converges to u 0 strongly in L p (Ω × Y m ). Let τ be a function in C ∞ 0 (Ω). We note that v η τ R → → v 0 τ (strongly)
. Thus, passing to the two-scale-cut-and-projection limit when η goes to zero in the product of v η τ and ψ m , we have that

(31) lim η→0 Ω ψ m x, Rx η v η (x)τ (x) dx = Ω×Y m ψ m (x, y) v 0 (x, y)τ (x) dxdy
We then pass to the limit when m goes to infinity (32)

lim m→∞ lim η→0 Ω ψ m x, Rx η v η (x)τ (x) dx = Ω×Y m u 0 (x, y) v 0 (x, y)τ (x) dxdy
where we have used that ψ m converges to u 0 strongly in L p (Ω × Y m ). Moreover, from the triangular inequality, we have (33

) Ω u η (x)v η (x)τ (x) dx - Ω×Y m u 0 (x, y) v 0 (x, y)τ (x) dxdy ≤ Ω u η (x) -ψ m x, Rx η v η (x)τ (x) dx + Ω ψ m x, Rx η v η (x)τ (x) dx - Ω×Y m u 0 (x, y) v 0 (x, y)τ (x) dxdy
Combining (32) and (33) we get

(34) lim sup η→0 Ω u η (x)v η (x)τ (x) dx - Ω×Y m u 0 (x, y) v 0 (x, y)τ (x) dxdy ≤ lim sup m→∞ lim sup η→0 Ω u η (x) -ψ m x, Rx η v η (x)τ (x) dx
It remains to prove that the right-hand side in (34) vanishes.

We first recall that τ ∈ C ∞ 0 (Ω) and we then invoke Hölder's inequality to deduce that ( 35)

Ω u η (x) -ψ m x, Rx η v η (x)τ (x) dx ≤ max x∈Ω | τ (x) | Ω u η (x) -ψ m x, Rx η v η (x) dx ≤ max x∈Ω | τ (x) | Ω u η (x) -ψ m x, Rx η p dx 1/p Ω | v η (x) | q dx 1/q ≤ C(Ω) u η -ψ m x, Rx η L p (Ω)
where in the last inequality we have used that v η is a bounded sequence in L q (Ω). We now invoke the Clarkson inequalities applied to functions u η and ψ m in L p (Ω):

(36)

1 2 p u η -ψ m x, Rx η p L p (Ω) ≤ 1 2 u η p L p (Ω) + ψ m x, Rx η p L p (Ω) - u η + ψ m x, Rx η 2 p L p (Ω)
for 1 < p < 2 and (37)

1 2 q u η -ψ m x, Rx η q L p (Ω) ≤ 1 2 u η p L p (Ω) + 1 2 ψ m x, Rx η p L p (Ω) q p - u η + ψ m x, Rx η 2 q L p (Ω)
for p ≥ 2

Passing to the 2-scale cut-and-projection limit in (36) (38)

lim sup η→0 u η -ψ m x, Rx η p L p (Ω) ≤ 2 p-1 u 0 (x, y) p L p (Ω×Y m ) + ψ m (x, y) p L p (Ω×Y m ) -2 p u 0 (x, y) + ψ m (x, y) 2 p L p (Ω×Y m )
for 1 < p < 2 and similarly for (37)

(39) lim sup η→0 u η -ψ m x, Rx η q L p (Ω) ≤ 2 u 0 (x, y) p L p (Ω×Y m ) + ψ m (x, y) p L p (Ω×Y m ) 1 p-1 -2 q u 0 (x, y) + ψ m (x, y) 2 q L p (Ω×Y m )
for p ≥ 2 where we have used also that 1/p + 1/q = 1 thus q/p = 1/(p -1).

We now note that u 0 (x, y) + ψ m (x, y) 2

L p (Ω×Y m )
→ u 0 (x, y) L p (Ω×Y m ) when m tends to infinity. Thus, taking the lim sup on m in both sides of (38) and (39), we are ensured that (40) lim sup

m→∞ lim sup η→0 u η -ψ m x, Rx η L p (Ω) = 0 for 1 < p < ∞
We have thus proved that the RHS of (34) vanishes combining ( 35) and ( 40).

(iii) We finally want to show (25). Assuming that u 0 ∈ L p (Ω, C (Y m )), we can replace ψ m by u 0 in the Clarkson inequalities (36) and (37), which leads to (41) lim sup

m→∞ lim sup η→0 u η -u 0 x, Rx η L p (Ω) = 0 for 1 < p < ∞ which is (25).
Remark 1.2. Our proof of proposition 1.3 follows closely that of Theorems 10 and 11 in [START_REF] Lukkassen | Two-scale convergence[END_REF].

We stress that some regularity is needed for u 0 in (iii) of proposition 1.3. We refer to Theorem 12 in [START_REF] Lukkassen | Two-scale convergence[END_REF] for a related result on regularity of two-scale limit that remains valid for two-scale-cut-and-projection limit (its proof is a straightforward extension of that in [START_REF] Lukkassen | Two-scale convergence[END_REF], but is rather technical and lengthy).

Remark 1.3. We point out that Proposition 1.3 i) does not hold if weak twoscale convergence u η k

R u 0 is replaced by distributional two-scale convergence (see definition 1.1). Indeed, the choice of space L q (Ω, C (Y m )) for test functions ϕ in (1.1) is essential. This is exemplified by the counter-example in [START_REF] Lukkassen | Two-scale convergence[END_REF] of a sequence (u η ) in L p (0, 1) defined by u η (x) = 1/η if 0 < x < η and u η (x) = 0 if η < x < 1. One can see that (6) is satisfied for ϕ ∈ D((0, 1); C ∞ (Y m )) and a two-scale limit u 0 (x, y) = 0. However, considering the test function g(x) = 1 which is in L q (0, 1), we get lim η→0 1 0 u η (x)g(x)dx = 1, and so (u η ) does not converge to u 0 (x, y) = 0 weakly in L p (0, 1).

Classes of functions such that u 0 x, Rx η

L p (Ω) → u 0 (x, y) L p (Ω×Y m ) are
said to be admissible for the two-scale (cut-and-projection) convergence. In particular, classes of functions in

L p (Ω, C (Y m )) (dense subset in L p (Ω × Y m )) are admissible.
In order to homogenize nonlinear PDEs with a monotone partial differential operator as in (3), we need to identify the differential relationship between χ and u 0 , given a bounded sequence (u η ) in W 1,p (Ω) (such that u η R u 0 and ∇u η R χ).

This problem was solved by Allaire in the case of periodic functions [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] and extended by Bouchitté et al. for quasiperiodic functions [START_REF] Bouchitté | Homogenization of dielectric photonic quasi crystals[END_REF] in W 1,2 (Ω) and revisited in [START_REF] Wellander | Two-scale cut-and-projection convergence; homogenization of quasiperiodic structures[END_REF][START_REF] Wellander | Homogenization of quasiperiodic structures and two-scale cut-and-projection convergence[END_REF].

Function spaces for cut-and-projection partial differential operators

To carry out the homogenization analysis of nonlinear PDEs defined on quasiperiodic domains, we need to pass to the limit when η goes to zero in gradient and divergence operators acting on solutions of PDEs. To do this we introduce some suitable function spaces and for this we will define differential operators acting on the R n -subspace in R m . Defined as in [START_REF] Wellander | Two-scale cut-and-projection convergence; homogenization of quasiperiodic structures[END_REF] they are given as

∇ R u(y) = grad R u(y) = R T grad y u(y) = R T ∇ y u(y) div R u(y) = R T ∇ y • u(y)
We define the following functions spaces associated with the differential operators defined above

(42) W p (grad R , Y m ) = u ∈ L p (Y m ) | grad R u ∈ L p (Y m ; R n ) (43) W p (div R , Y m ) = u ∈ L p (Y m ; R n ) | div R u ∈ L p (Y m ) and (44) W p (div, Ω) = u ∈ L p (Ω; R n ) | div u ∈ L p (Ω
We have the following integration by parts type generalization to the L p case of Lemma 6 given in the L 2 setting in [START_REF] Wellander | Homogenization of quasiperiodic structures and two-scale cut-and-projection convergence[END_REF] Lemma 2.1 (Green's Identity). It holds that

(45) - Y m RR T ∇ y • φ(y) θ(y) dy = Y m φ(y) • RR T ∇ y θ(y) dy for φ ∈ W q (div R , Y m ) and θ ∈ W p (grad R , Y m ), 1/p + 1/q = 1.
Proof The proof relies on standard matrix operations and the well known extension from W 1,2 Sobolev spaces to the setting of W 1,p and W 1,q -duality pairing [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]. The periodic boundary conditions imply (46)

- Y m RR T ∇ y • φ(y) θ(y) dy = - Y m R T ∇ y • R T φ(y) θ(y) dy = - Y m ∇ y • RR T φ(y) θ(y) dy = Y m RR T φ • ∇ y θ(y) dy = Y m φ(y) • RR T ∇ y θ(y) dy for any pair of functions for φ ∈ W q (div R , Y m ) and θ ∈ W p (grad R , Y m ), 1/p + 1/q = 1.
This motivates us to make the following decomposition as in the L 2 -case in [START_REF] Wellander | Homogenization of quasiperiodic structures and two-scale cut-and-projection convergence[END_REF]. We decompose W 1,p (Y m ) into two spaces, (47)

W 1,p (Y m ) = X p ⊕ X ⊥ p where (48) X ⊥ p = u ∈ W 1,p (Y m )|RR T ∇ y u = 0 and (49) X p = u ∈ W 1,p (Y m )| I m -RR T ∇ y u = 0 Remark 2.1. Note that the projection of a vector v in R m on R n , R T v = 0,
where 0 is the zero vector in R n , implies that v is orthogonal to the hyperplane, i.e., orthogonal to R n . It follows that RR T v = 0, where 0 is the zero vector in R m . Vectors w in R m , orthogonal to v, satisfies I m -RR T w = 0. We conclude that X p contains all functions in W 1,p (Y m ) whose gradients have all their components in the hyperplane, R n , which means that X p can be identified with W p (grad R , Y m ).

Compactness, homogenization and corrector results

Proposition 3.1. Let {u η } be a uniformly bounded sequence in W 1,p (Ω). Then there exist a subsequence {u η k } and functions u ∈ W 1,p (Ω) and

u 1 (x, y) ∈ L p (Ω, X p ) such that (50) u η k R → → u(x), grad u η k R grad u(x) + grad R u 1 (x, y), η k → 0 Remark 3.1. Note that u 1 (x, y) ∈ L p (Ω, X p ) implies that grad R u 1 (x, y) ∈ L p (Ω, L p (Y m ; R n ))
and that we don't have enough regularity to ensure that u 1 belongs to L p (Ω, W 1,p (Y )). Indeed, we cannot say anything about the regularity of u 1 in the direction orthogonal to the hyperplane. However, in the decomposition in (47), the gradient of the potential in the "direction " of the hyperplane can be obtained from the gradient of the potential via a rotation of coordinate system. Further note that unlike in [START_REF] Wellander | Two-scale cut-and-projection convergence; homogenization of quasiperiodic structures[END_REF][START_REF] Wellander | Homogenization of quasiperiodic structures and two-scale cut-and-projection convergence[END_REF], in the proof below we use the notion of 2-scale cut-and-projection convergence in distributional sense (Definition 1.1), and not in weak sense (Definition 1.2).

Proof The first assertion follows by the compact embedding of L p (Ω) in W 1,p (Ω), Propositions 1.1 and Definition 1.3. Note that R∇u η (x),

x ∈ Ω is uniformly bounded in L p (Ω; R m ). Let ϕ ∈ D(Ω; C ∞ (Y m )) m . We have the following identities (51) lim η→0 Ω R∇u η (x) • ϕ x, Rx η dx = lim η→0 Ω ∇u η (x) • R T ϕ x, Rx η dx = -lim η→0 Ω u η (x) ∇ x • R T ϕ x, Rx η + η -1 R T ∇ y • R T ϕ x, Rx η dx and (52) R T ∇ y • R T ϕ x, Rx η = RR T ∇ y • ϕ x, Rx η
Multiplying both sides in (51) with η and Lemma 2.1 gives the limit (53)

0 = Ω×Y m u 0 (x, y) R T ∇ y • R T ϕ (x, y) dxdy = Ω×Y m u 0 (x, y) RR T ∇ y • ϕ (x, y) dxdy = Ω×Y m RR T ∇ y u 0 (x, y) • ϕ (x, y) dxdy for all ϕ ∈ D(Ω; C ∞ (Y m )) m .
The interpretation of ( 53) is that u 0 ∈ X ⊥ p , i.e., the gradient ∇ y u 0 (x, y) has no component in the hyper plane in R m defined by R : R n → R m . Indeed, we conclude that the potential u 0 (x, y) = u(x), is a function of x only due to the compact embedding of L p (Ω) in W 1,p (Ω) and that the two-scale limit equals the strong limit, if it exists. Next, let ϕ ∈ D(Ω; C ∞ (Y m )) m , and ψ ∈ D(Ω; C ∞ (Y m )) n . We have three limits We deduce, due to orthogonality in the dual pairing sense (45), that there exists u 1 ∈ L p (Ω; X p ) such that χ0 (x, y) = R∇u(x) + RR T ∇ y u 1 (x, y)

lim η→0 Ω R∇u η (x) • ϕ x, Rx η dx = -lim η→0 Ω u η (x) ∇ x • R T ϕ x, Rx η dx = - Ω×Y m u(x) ∇ x • R T ϕ (x, y) dxdy = Ω×Y m R∇u(x) • ϕ (x,
where u(x) = Y m u 0 (x, y)dy, where u 0 (x, y) is the two-scale cut-and-project limit with respect to R. Next, we have the corresponding two-scale cut-and-project limit of the divergence (64)

lim η→0 Ω ∇ • u η (x)ψ x, Rx η dx = Ω×Y m χ 0 (x, y)ψ(x, y) dxdy , ∀ψ ∈ L q (Ω, C (Y m ))
It follows, after an integration by parts (twice) that (65)

0 = lim η→0 η Ω ∇ • u η (x)ψ x, Rx η dx = Ω×Y m div R u 0 (x, y)ψ(x, y) dxdy , ∀ψ ∈ L q (Ω, C (Y m ))
which proves (62). Define a function as the difference of the two-scale and the weak limits, i.e., f (x, y) := χ 0 (x, y) -∇ • u(x). We have f (x, •) ∈ L p (Y m ). Lemma 3.1 implies that there is a unique grad R θ(x, •) in L q (Y m ; R n ) that solves (60), i.e., θ(x, •) ∈ X q , defined in (49). Next, define u 1 (x, y) := -a (grad R θ(x, y)) with a as in Lemma 3.1. We get

χ 0 (x, y) = div u(x) + f (x, y) = div u(x) + div R u 1 (x, y) ∈ L p (Ω × Y m )
which completes the proof.

Homogenization of a quasiperiodic heterogeneous nonlinear electrostatic problem

Let us now consider the quasiperiodic heterogeneous nonlinear electrostatic problem (3). Standard estimates yield solutions that are uniformly bounded in W 1,p 0 (Ω) with respect to η. We can now state the main homogenization result.

Theorem 4.1. Let {u η } be a sequence of solutions to (3). The whole sequence {u η } converges weakly in W 1,p 0 (Ω) to the solution {u} of the homogenized equation

(66)    -div Y m σ x, y, ∇u(x) + R T ∇ y u 1 (x, y) dy = f (x) , x ∈ Ω u| ∂Ω = 0 where R T ∇ y u 1 (u 1 ∈ X p ) is the unique solution of the local equation (67) -div R σ x, y, ∇u(x) + R T ∇ y u 1 (x, y) = 0 x ∈ Ω
Proof From the a priori estimates of sequences u η and σ

η := σ x, Rx η , ∇u η (x) , there is subsequence such that u η k R → → u(x), grad u η k R grad u(x)+grad R u 1 (x, y)
and σ η R σ 0 (x, y), when η k → 0. Since f + divσ η = 0, equation (62) yields div R σ 0 (x, y) = 0 and f (x) + div x Y m σ 0 (x, y)dy = 0. We now need to obtain an explicit expression for σ 0 (x, y) in terms of σ, u and u 1 . Following, e.g., [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] we introduce a test function ψ η (x) = ∇{u(x) + ηφ 1 (x, Rx η )} + tφ(x, Rx η ) where t > 0, φ and φ 1 are admissible test functions. This ensures that ψ η k R grad u(x) + grad R φ 1 (x, y) + tφ(x, y). Since σ is strictly monotone, we have ( 68)

Ω σ η -σ x, Rx η , ψ η • (∇u η -ψ η ) dx ≥ 0 i.e., (69)
Ω -divσ η u η -σ x, Rx η , ψ η • ∇u η -σ η ψ η + σ x, Rx η , ψ η ψ η dx ≥ 0
Using (3), passing to the two-scale limit and using the strong limit to get u yields

Ω Y m {f (x)u(x) -σ (x, y, ψ 0 (x, y)) • (grad u(x) + grad R u 1 (x, y)) -σ 0 (x, y)ψ 0 (x, y) + σ (x, y, ψ 0 (x, y)) ψ 0 (x, y)} dxdy ≥ 0
This equals, after a few integration by parts,

Ω Y m {f (x)u(x) + div σ (x, y, ψ 0 (x, y)) u(x) -σ (x, y, ψ 0 (x, y)) • grad R u 1 (x, y) + div σ 0 (x, y)u(x) + div R σ 0 (x, y)φ 1 (x, y) -σ 0 (x, y)tφ(x, y)--div σ (x, y, ψ 0 (x, y)) u(x) -div R σ (x, y, ψ 0 (x, y)) φ 1 (x, y)+ σ (x, y, ψ 0 (x, y)) tφ(x, y)} dxdy ≥ 0

The first terms in the first two rows cancel each other due to the statements above. We also note that the middle term in the second row vanishes. The second term in the first row is canceled by the first term in the third row. Taking a sequence of functions grad R φ 1 that converges strongly to grad R u 1 in L p (Ω, L p (Y n ) n )), yields

Ω Y m {σ(x, y, grad u(x) + grad R u 1 (x, y) + tφ(x, y)) • grad R u 1 (x, y) + -σ 0 (x, y)tφ(x, y) -σ(x, y, grad u(x) + grad R u 1 (x, y) + tφ(x, y)) • grad R u 1 (x, y)+ σ(x, y, grad u(x) + grad R u 1 (x, y) + tφ(x, y))tφ(x, y)} dxdy ≥ 0

The first row cancels the second term in the second row. We divide the two terms left by t > 0 and send t to zero and obtain If follows from the fact that ψ η are admissible test functions that the left hand side of (73) goes to zero in L p (Ω; R n ).

Concluding remarks

We have applied two-scale cut-and-projection convergence to a canonical nonlinear electrostatic problem for quasiperiodic structures generated by a periodic geometry in a higher dimensional space. Compared with earlier work on homogenization of almost periodic monotone operators [START_REF] Braides | Correctors for the homogenization of almost periodic monotone operators[END_REF][START_REF] Nguetseng | Sigma-convergence of nonlinear parabolic operators[END_REF], our annex problem has a simpler, less abstract structure, and should therefore facilitate its numerical implementation in a variety of problems of physical interest, such as in electromagnetism [START_REF] Wellander | Homogenization of some linear and nonlinear partial differential equations[END_REF], where intriguing features have been observed, such as transmitted femtosecond pulses developed a trailing diffusive exponential tail that led to some controversy [START_REF] Ledermann | Multiple scattering of light in three-dimensional photonic quasicrystals[END_REF]. We further note that our study can be adapted to the nonlinear elasticity case [START_REF] Castaneda | The overall constitutive behaviour of nonlinearly elastic composites[END_REF], whereby σ would denote a rank-4 elasticity tensor.

  R∇u η (x) • ϕ x, Rx η dx = Ω×Y m χ0 (x, y) • ϕ(x, y) dxdy (55) lim η→0 Ω R∇u η (x) • ϕ x, Rx η dx = lim η→0 Ω ∇u η (x) • R T ϕ x, Rx η dx = Ω×Y m χ 0 (x, y) • R T ϕ(x, y) dxdy = Ω×Y m Rχ 0 (x, y) • ϕ(x, y)Ω×Y m χ0 (x, y) • ψ(x, y) dxdy We find that χ0 (x, y) = Rχ 0 (x, y) and χ0 (x, y) = R T χ0 (x, y) = χ 0 (x, y). Next, choosing test functions ϕ ∈ D(Ω; C ∞ (Y m )) m , such that RR T ∇ y • ϕ(x, y) = 0 in (51) with (52) gives (57)

  y) dxdyHence due to (54) and (57) we have forϕ ∈ D(Ω; C ∞ (Y m )) m , such that RR T ∇ y • ϕx, y) -R∇u(x)) • ϕ (x, y) dxdy = 0

0 Proof

 0 x, y, grad u(x) + grad R u 1 (x, y)) -σ 0 (x, y)] φ(x, y)dxdy ≥ 0 for all admissible test function, e.g., φ ∈ D(Ω; C ∞ (Y m )). It follows that σ 0 (x, y) = σ(x, y, grad u(x)+grad R u 1 (x, y)). Uniqueness of the solution of the limit equation (see e.g.[START_REF] Lions | Quelques méthodes de résolution des problemes aux limites non linéaires[END_REF][START_REF] Wellander | Homogenization of some linear and nonlinear partial differential equations[END_REF]) implies that the whole sequence converges. Proposition 4.1 (Correctors). If we assume that u 1 (x, y) is smooth and σ is uniformly monotone, thenlim η→0 ∇u η (x) -∇ u(x) + ηu 1 x, Rx η L p (Ω;R n ) = Considering ψ η (x) = ∇ u (x) + ηu 1 x, Rx ηand using that σ is uniformly monotone yields (73)Ω σ η -σ x, Rx η , ψ η • (∇u η -ψ η ) dx ≥ c Ω | ∇u η -ψ η | p dx

Acknowledgement

The work of E.C. was supported by NSF grant DMS-1715680. S.G. acknowledges funding from the UK Engineering and Physical Sciences Research Council (Grant No. EP/T002654/1)

We conclude that the limit of the gradient in (56) becomes (59) χ0 (x, y = R T χ0 (x, y) = R T R∇u(x) + RR T ∇ y u 1 (x, y) = ∇u(x) + R T ∇ y u 1 (x, y) = grad u(x) + grad R u 1 (x, y) which completes the proof.

We define a strictly monotone operator a, which satisfies the following assumptions, i)-iv

We will have use of the following Lemma when characterizing the two-scale limit of divergences. Lemma 3.1. Let 1 < p < ∞, 1/p + 1/q = 1 and assume the operator a satisfies assumptions i)-iv) above and that f (x, •) ∈ L q (Y m ). The equation

Proof The proof follows from Browder (1963) and Minty (1963), e.g., see [START_REF] Lukkassen | Two-scale convergence[END_REF], page 62. Proposition 3.2. Let {u η } be a uniformly bounded sequence in W(div, Ω). Then there exist a subsequence {u η k } and functions in u 0 ∈ W p (div, Ω; L p (Y m )) and

Proof The proof follows the lines of Lemma 5 and Proposition 6 in [WGC19] with appropriately changed function spaces. Let φ ∈ L q (Ω) and ψ ∈ L q (Ω, C (Y m )). We have the weak limit of the divergence, (63)