Change of boundary: towards a mathematical foundation of global gravity models

Change of boundary is a method that iteratively downward continues data from a star-shaped boundary to a regular surface, such as the sphere or the ellipsoid of revolution, and then solves a boundary value problem using the downward continued data on the changed boundary. Although the method belongs to concepts involved in the development of recent high-degree Earth’s gravity models, it is still unclear whether the iterations converge to the true solution for boundaries as complex as, for instance, the Earth’s surface. In this paper, we revisit the method and show that convergence in terms of the Cesaro limit can be achieved under the assumption that the operator performing the iterations is non-expansive. The validity of the assumption is, however, still not proved. Therefore, we examine the hypothesis numerically using boundaries of various complexity. We start with a simple synthetic topography defined by a Legendre polynomial and move to more realistic finite-degree shapes of the asteroids Bennu and Eros. The numerical experiments indicate that the assumption is valid as long as the boundary deviates not too far from a sphere and the truncation degree of the gravity model is not too high (the experiments with the synthetic topography and Bennu). Otherwise, the hypothesis seems to be false (the Eros case). Finding an analytical condition to separate between shapes for which the change of boundary method converges/diverges remains an open issue.


Introduction and notation
The determination and use of global gravity models has been the story of a success of physical geodesy in the last decades. A global model of the external anomalous potential of the Earth T (r , σ ), in this context, is just a harmonic function represented by a truncated series of solid spherical harmonics (Pavlis 2013;Sansò 2013b) where r is the spherical radius and is a pair of spherical angles, colatitude and longitude. Y nm (σ ) are fully-normalized real spherical harmonics and u nm are spherical harmonic coefficients. N is the maximum (or truncation) degree of the model. R is in principle any radius, essentially defining the length unit of the problem; in fact, it is obvious that, under the transformation u(r , σ ) is left invariant.
In the literature, we find the choice R = a, the equatorial radius of the ellipsoid, or R = mean radius of the Earth. In this paper, we will use for R any radius of a sphere S R , internal and separated from the Earth surface S ≡ {r = R σ }, that we assume to be a star-shaped and Lipschitz surface. Let us notice here too, that we will use the notation for the domain exterior to S. Furthermore, the notation R σ , which is a shorthand for a function of σ , should not be confused with R, which is the constant radius of a sphere.
If we put (i.e. the sup of the Bjerhammar radius and inf of the Brillouin radius), we will assume that The reason why R is so loosely defined is that, with the up and down transformation (3), we leave u(r , σ ) invariant and what really matters is that the global model (1) needs to be a "good" approximation of the actual anomalous potential. On this point one has to be very clear. What we are looking for is a theory that certifies that our global potential (1) is a good approximation of the anomalous potential T (r , σ ) in (or in some cases on S, too), in the sense that in this domain lim N →∞ u(r , σ ) = T (r , σ ) in some topology for a space of functions harmonic in .
In this paper, we will use exclusively the harmonic Hardy space HL 2 with norm defined by Let us note that we find it convenient to use in the paper a normalized measure on the unit sphere, that we will denote dσ = 1 4π sin ϑ dϑ dλ.
Fatou's lemma and Cimmino's theory (Axler et al. 1992 revisited in Sideris 2017 andVenuti 2001) guarantee that the limit (7) exists and in fact HL 2 is in a unitary correspondence with L 2 (σ ), namely the Hilbert space of functions of σ with norm of course on condition that (7) be bounded. The choice of this topology is because it is strongly set in geodetic tradition and, when discretized, it leads to least squares approximation, which is a universal tool in the arsenal of geodetic methods. In addition, the theory of (linearized) geodetic boundary value problems has shown that once a datum of free air gravity anomaly or disturbance is given on S, belonging to L 2 (σ ), then a solution (or anomalous potential T ) exists, is unique and continuously dependent on the data in L 2 (σ ) norm (Sansò 2013a). This fact together with a theorem of completeness of the base functions where which is a particular case of Runge-Krarup theorem (Krarup 1975;Sansò 2013c), says that the problem of the approximation of T (r , σ ) by a global potential of the form (1) is well-posed in L 2 (σ ). Yet the way in which the most important global models at high resolution (N 2159) are actually computed (see for instance Pavlis 2013) is not by the use of least squares directly applied to the original data, but rather the data are moved ("downward continued") from the surface down to the ellipsoid, which is then perturbatively approximated by a sphere, and finally the global model is computed as the solution of the boundary value problem from this easy to treat "changed boundary".
Residuals are then computed at the level of S, where data are given, and the process is iterated. Indeed the real way in which a global model is computed includes a number of very important intermediate steps, including the transformation of marine altimetric data into gravity anomalies and the assimilation of satellite only gravity models (see Pavlis 2013;Rapp 1993).
In particular, this last step has a large influence in the estimation of harmonic coefficients up to some intermediate degree (a few hundreds) and has a fundamental role in the correction of biases present in the ground data, due to the nonhomogeneous historical measurement process (height datum problem; Sansò and Venuti 2001;Rapp 1993).
Yet for the high-frequency (degree) coefficients the determination process is essentially outlined above. Therefore, the question arises on how good is the result obtained. No doubt the answer has to be given by testing the estimated field on independent data or by cross-validation.
This has been quite successfully done (e.g., Pavlis et al 2012) and therefore it is clear, on a scientific ground, that a mathematical proof of convergence of the iterations is needed to explain/understand such a success. This is how the idea of the change of boundary (CoB) theory was born. Up to now the results obtained are not complete, because of a hard obstacle in proving that a certain operator is non-expansive. Some steps forward have been done by conjecturing that this assumption is correct, based on a perturbative argument (Sansò and Sideris 2017); yet a rigorous proof is still lacking.
In such a situation it was deemed useful by the authors to set up controlled simulations to test at least numerically the validity of the conjecture. This is what is done in this paper, together with some mathematical refinements, showing the results that can be obtained if the mentioned assumption is correct.
In particular a verifiable sufficient condition is given, to guarantee the convergence if we pose an a priori hypothesis that in real world, or better in our model world, both the anomalous potential and the shape of S are expressed by spherical developments with known maximum degree. Therefore, we want to stress again that the purpose of the paper is not to propose CoB as a new method to estimate global models, but rather to set up a formal frame in which to understand the success of existing methods.
Concluding this section, it is worth mentioning that the item of convergence of global models has been tackled by several authors, also outside the geodetic literature (e.g., Takahashi and Scheeres 2014;Hu and Jekeli 2015;Reimond and Baur 2016;Sebera et al. 2016;Hirt and Kuhn 2017). However, the purpose of this paper is slightly different in that we try to assess the extent at which our present theoretical understanding of the problem is reliable; in particular this will be done by examining the norm of a certain matrix operator, rather than following the sequence of models of increasing degree.

Foundation of CoB
In order to set up the formal definition of CoB we need to recall first of all that the linearized geodetic boundary value problem amounts to determining T (r , σ ) in given, depending on what we consider as known on the boundary, on S ≡ {r = R σ }; here h is as usual the ellipsoidal height, γ the modulus of the normal gravity vector, δg the gravity disturbance and g the gravity anomaly. Let us recall that the CoB has been first defined in Sansò (1993) and its theory more carefully analysed in Sansò and Sideris (2017). The spherical approximation of (12) (see Sansò and Sideris 2017) is just the so-called simple Molodensky problem and its solution is the first step for the solution of (12). The point is that in either cases (13) the boundary value problem can be reduced to the solution of a Dirichlet problem; in fact, one has that have to be harmonic too (Heiskanen and Moritz 1967;Sansò and Sideris 2017) in and they satisfy one or the other boundary condition For the sake of simplicity let us choose one of the two, for instance the first that has no singularity at degree 1, as it happens to the second. So the problem of finding T harmonic in ≡ {r > R σ }, satisfying the boundary condition (15), is equivalent to searching for the solution u(r , σ ) of the Dirichlet problem with boundary values (15), and then retrieving T (r , σ ) by integration of −u(r , σ )/r along the radius (Sansò 2013a). The geodetic practice is then to "lower" f (σ ) from S to S R ; this step can be done in different ways, for instance by collocation or simply by a truncated Taylor series (Martinec 1998;Pavlis 2013). Whatever it is, we call PB (pull back) the operator that lowers f (σ ) to S R . After that the Dirichlet problem is solved for the pull back PB f by Poisson formula and the residuals between f and the upward continuation of PB f are computed. Would the residuals be zero, we could say that the problem has been solved. Since this is not the case, we will see that it cannot be, the same procedure is applied to the residuals and the same approach is repeated iteratively.
Then all contributions from the iterations are added to find hopefully a series convergent, in some sense, to the correct solution. The situation is formalized in the scheme of Fig. 1 where PB represents the pull back operator, denotes the Poisson operator, harmonically propagating the boundary values on the sphere S R to the outer domain R (i.e. the exterior of the sphere S R ) and is the trace on S operator, namely given that S ≡ {r = R σ }, is the operator that, for any continuous u(r , σ ) defined in R , computes Furthermore, we call U the upward continuation operator or, more precisely (Heiskanen and Moritz 1967) ( σ σ = [R 2 + R 2 σ − 2 R R σ cos ψ σ σ ] 1/2 , ψ σ σ = spherical angle between σ and σ ).
Notice that the Poisson formula (19) is written without the factor (4π) −1 which is included by definition in the area element dσ (see 8). Since S is star-shaped by hypothesis, U transforms functions of σ (i.e. defined on the unit sphere) into functions again defined on the unit sphere. Let us further observe that, if S ⊂ R and as we assume, any function harmonic in R is continuous there, so that (19) is well-defined. Moreover, in Fig. 1 we use D to denote the Dirichlet solver that, given a function f (σ ) on S, provides the function w(r , σ ), harmonic in and agreeing with f (σ ) on S. The scheme of Fig. 1 implies the following relations Combining the three relations one gets: So, putting one can write implying, also recalling (25), that On the same time, iterating (22) we see that Therefore, a first conclusion of our construction is that if in some topology we can claim that then we have found a sequence F (σ ), defined only by applications of PB and U (i.e. Poisson formula), that in the same sense (topology) tends to f (σ ).
In addition, still looking at Fig. 1, one realizes that which are the solutions of Dirichlet problem agreeing with f (σ ) on S, have to be just the restriction of u −1 (r , σ ) to ; in fact, these are harmonic in R > and they are equal to f (σ ) on S, by definition of f (σ ). Therefore, if the topologies used are such that we can conclude too that namely we have found the approximate solution we are looking for.
By expressing u k (r , σ ) in terms of solid spherical harmonics, what is natural since they are given by Poisson formula, we will also find convergence theorems for global models.
Remark 1 It is important to underline again that the operators PB and U = (see 19), which are the constituents of the iterative process, can be interpreted as transforming functions of σ only into functions of σ again so they can be thought of as transformations between some functional space into itself. In this paper this space, i.e. the topology in which we embed the problem, is simply L 2 (σ ), i.e. the Hilbert space with norm (31), with W (r , σ ) → w(r , σ ) interpreted as pointwise convergence or uniform convergence on compact subsets of , holds, in force of Runge-Krarup's theorem in the simplified form as presented on p. 137 of Sansò (2013b) (see also Krarup 1975, [p. 254]).
Before developing our analysis, we want to recall a first choice done for the operator PB: We justify such a position by examining the two sphere example, first presented in Sansò (1993). We are going to see in fact that any choice of PB as a Taylor formula in radial direction has essentially not very different characteristics than (35) and therefore this is the simplest assumption. Nevertheless, as we will see in Sect. 6, the alternative has been explored too, in an attempt to improve our results. Yet we will find that even more complex gravity field models produce results similar to those obtained by next example.
Since it will be used later on, we give also a spectral representation of PB from (36), when acting on a potential harmonic down to S R , namely Example 1 (Two spheres S R , S R 0 ) In this example we examine the case that S = S R 0 , a sphere of radius R 0 > R. This is because, despite the elementary character of this case, we can work it out analytically and read directly how the CoB scheme works. We put let us notice that in this example the three quantities q 0 , h, η are constant. We test the CoB iterations with the two choices In both cases we can write the spectral representation of PB as Moreover, the spectral representation of U is just and so with c n one of the two choices in (40). Clearly if k n < 1, the operator K as well as all K for any are continuous operators in L 2 (σ ).
Since q 0 < 1, it is obvious that k n < 1 when c n = 1, so we have to look into (42) when c n = 1 + h R 0 (n + 1). Considering that we easily get, for the second case in (40), Summarising, we have for k n the values corresponding to the first of (40), and corresponding to the second. As we see, we have It is not difficult to verify that the same would happen even if we pushed the Taylor expansion in PB to a higher order. Since the convergence of the iterative scheme critically depends on the behaviour of k n , which are the eigenvalues of K corresponding to the eigenfunctions Y nm (σ ), and the relation (47) holds, it is enough here to consider the case of k n (45) in order to prove convergence of the CoB scheme.
Two important comments are in order.

Remark 2
We have to underline that the forms (45), (46) of the eigenvalues of U , together with relations (47), (48) show that because 1, the spectral radius of K , is also equal to its operator norm in L 2 (σ ). This will be a hard point when we will try to extend this analysis to a general surface S different from a sphere S R 0 . On the other hand, if we truncated the operator K to some maximum degree N , then we would have and the convergence will immediately descend from the fact that K becomes a contraction. Also this remark will be used in a more general setting.

Remark 3
Considering Example 1, we have to understand that the coefficients { f nm } of f (σ ) refer to a function, the datum, which is given on the surface S R 0 . So the convergence of s to 0 leading to the convergence of the method is for functions defined on S R 0 and W (r , σ ) tend to w(r , σ ) in On the other hand, if we call i.e. the sum of the residuals seen at the level of S R (recall that R < R 0 ), we see from (23) that Now the eigenvalues of U are c n = q n+1 0 , so they are positive tending to 0 for n → ∞. It follows that U is a compact operator and then the convergence of F (σ ) to f (σ ) in L 2 (σ ) does not imply that V (σ ) is convergent too to some function in L 2 (σ ). Would this be true, we could always downward continue f (σ ) from S R 0 to S R , which is known to be false.
A final word of caution, before closing the section, is that although the choice PB = I makes the problem more comfortable to analyse, in practice a radial Taylor expansion should always be applied before going to numbers (see, e.g. Sansò and Sideris 2017;Martinec 1998). correct or not. Whatever we choose as PB, we have so it is important to get a quite clear understanding of the properties of U . They are summarized below in a list from a) to e). Since U can be represented as an integral operator, as in (19), most properties of U are just well-known classical properties of such operators, so they will be only formulated with reference to literature. We would like to underline that such a matter is well-known in mathematics since the times of Fredholm, more than one century ago, and is presented in classical old books of Analysis, like Goursat (1964).
(a) U is a compact operator. According to (19), U is an integral operator with kernel It is enough to use the identity to realize that Similarly, one has so that the following (very rough) upper and lower bounds hold for U (σ, σ ): So U (σ, σ ) is a bounded strictly positive function. Clearly U (σ, σ ) is also continuous in its arguments, in fact Lipschitz continuous, as it is R σ by hypothesis. Whence, since the unit sphere is a compact set with measure 1, by the hypothesis done on dσ , we have a fortiori The first conclusion then is that U is a continuous operator L 2 (σ ) → L 2 (σ ). Moreover, U is (sequentially) completely continuous, namely it transforms any weakly L 2 (σ )-convergent series into a strongly convergent one.
where q σ is defined by (11) and P n are the ordinary Legendre polynomials. In fact, (64) is nothing but the well-known series the inclusion being dense and, as it is proved in Theo- . This is as a matter of fact a by-product of Runge-Krarup's theorem. (c) The compact operator U has a representation in terms of singular elements where This representation is also standard (see Appendix of Sansò and Sideris 2017). Since the kernel of the operator it is clear that the kernel U (σ, σ ) is given by To verify that the series (73) is not only formal one can observe that indeed which shows that the series is L 2 (σ, σ ) convergent. Further on, by using the orthonormality of {ξ i }, we can write, ∀u ∈ L 2 (σ ) Since also {ϕ i } is orthonormal in L 2 (σ ), one has Therefore, if we put λ i in decreasing order, one can state that Finally, we observe that since by (68) with dense embedding, and since R(U ) is densely contained in L 2 (σ ) by point b), we have that also span{ξ i } is dense in L 2 (σ ). Since {ξ i } is orthonormal we conclude that it is an orthonormal basis of L 2 (σ ).
(d) U is an invertible operator, namely the equation v = U u = 0 ( 8 0 ) has only the solution u = 0. In fact, consider the potential which is harmonic in R . By (80) v(R σ , σ ) ≡ 0, so that v ≡ 0 in . But then, by the unique continuation property, v ≡ 0 in R , too; so as it was to be proved. We notice that, from (71), no λ i can be equal to zero; i.e. λ i is an infinite bounded sequence tending to zero through positive values. Moreover, since and since by the above N (U ) = ∅, we conclude that therefore also {ϕ i } is an orthonormal basis of L 2 (σ ). (e) As a compact operator we know that U has only a point spectrum of eigenvalues {μ i } such that Since μ i can never be zero, {μ i } is an infinite sequence.
Since U (σ, σ ) is real and non-symmetric, it is clear that in general μ i are complex and they come in couples, i.e. if μ i is an eigenvalue, μ i is an eigenvalue too, unless they are both real. Moreover, we know that the spectral radius ρ of U is also the norm of U , i.e.
Since {μ i } is a bounded sequence with zero as unique accumulation point it has to be, by suitably choosing μ 1 , Therefore, there is an eigenfunction p 1 (σ ) such that Now (88) tells us that p 1 (σ ) has to be a continuous function, so that we find On the other hand, since > 0, as it is clear also from (64). Therefore (89) yields namely So under the geometrical hypotheses done, particularly (5), U is a contraction and K = I − U is an invertible operator.
Remark 4 Let us observe that the operator K = I −U , entering into our theory with the choice PB = I , has indeed eigenvalues (1 − μ i ) and its spectral radius is On the other hand, since μ i → 0 for i → ∞, we must have Therefore, considering that, |μ i | ≤ |μ 1 | = q 0 < 1, the only possibility to have K = 1 is that the complex values μ i belong to the filled area in Fig. 2. In particular, any real negative value of μ i would imply K > 1.

Remark 5
We underline that the abstract representation of U by (68) and of its kernel U (σ, σ ), (73), is once more an old achievement of Analysis, going back to Schmidt and Picard (Goursat 1964). The form (68) in particular provides an explicit expression for U −1 , i.e.
attention has to be paid to the fact that (72) cannot be represented by a kernel in L 2 (σ, σ ), because λ i → 0. In fact U −1 , as the inverse of a compact operator, is an unbounded has a solution that in terms of components can be written as showing that u ∈ L 2 (σ ) only if ( Goursat 1964, p. 158) So it is clear that not for every v ∈ L 2 (σ ) there is a solution u ∈ L 2 (σ ) of (96); only if (98) is verified, which is also known as Picard's condition, can u be "reasonably" computed. This is just another way to claim that the downward continuation operator is unbounded in L 2 (σ ) and therefore the downward continuation operation is improperly posed in L 2 (σ ).
seen as acting on the whole L 2 (σ ), the operator K = I − U cannot have a norm smaller than 1; therefore, the question is whether we have if on the contrary we would decide that K > 1, the whole CoB concept would be destroyed. Yet in this section we would like to analyse the consequences of (100). In other words, would K be equal to 1, what happens to the residuals s given by (28)? As already recalled, (100) is a "reasonable" hypothesis descending from the sandwich conjecture, namely that calling h nm the eigenvalues of one has therefore, would (102) be correct, we should have The relation (102) is only guessed on the basis of a perturbative argument (see Sansò and Sideris 2017, Appendix). We shortly recapitulate the reasoning for the right inequality in (102); a similar reasoning holds for the left inequality.
Proof Let S R be the sphere completely internal to S, so that R < R 0 < R (recall the definition 4) and S is totally inside the layer between S R 0 and S R . Let us put where as always q σ = R R σ , q 0 = R R 0 and q σ < q 0 . Then Indeed U 0 is selfadjoint with eigenvalues h 0nm = q n+1 0 and eigenvectors {Y nm (σ )}.
The authors want to comment once more that the above proof is only heuristic in character, and that it holds under the condition that (105) be satisfied, i.e. that the linearization of U with respect to U 0 gives a good approximation of this operator.
So we assume that (100) is true, namely that K is a nonexpansive operator. An immediate consequence is that s is a bounded sequence and therefore a subsequence of it is weakly convergent; yet the point is whether such a limit is 0 or not and whether 0 is the limit for the whole sequence. Even more we would prefer to have a stronger result, namely that s → 0 in the L 2 (σ ) norm. This is not possible in general, yet a close result can be obtained by using the concept of Cesaro limit.

Definition 1
We say that the sequence of real numbers {x n } has the Cesaro limit x if the Cesaro means In this case we write Remark 6 It is obvious that such a definition can be generalized to any topological space, where the definition of the limit of a sequence is specified by the topology. So, for instance, in a Hilbert space X we say that Similarly, endowing X of the weak topology, we say that Supplied with this definition, we turn our attention to the sequence {s } ∈ L 2 (σ ) and recalling (28) we see that their Cesaro means can be written as It is therefore convenient to define the mean operators it is obvious that, like K , all the mean operators are nonexpansive, namely Considering the sequence M N ( f ) ≡ C N ({s n }), we get an immediate result recurring to the ergodic theory, more specifically to the Mean Ergodic Theorem (Yosida 1995). i.e.
Proof We adapt here the very general proof given in Yosida (1995) to our much simpler case. We first notice that the following identity holds Now, given (120) we know that ∀ f ∈ L 2 (σ ) fixed, the sequence M N ( f ) is bounded by g 2 ; so there is a subsequence {N j } such that for some g ∈ L 2 (σ ). Therefore by (124) and on the same time (K is continuous) (127) so we must have because of property d) in Sect. 3. We conclude that Since the same reasoning can be repeated for any other subsequence {N k }, we can conclude that for the whole sequence Hence Now take any w ∈ R(U ), i.e and notice that Since M N (u) − → 0 and U is compact by property a) from Sect. 3, we find On the other hand we know that R(U ) is dense in L 2 (σ ) by property b) from 3; then given any f ∈ L 2 (σ ) and any ε > 0 we can find w ∈ R(U ) such that Then we have Therefore, by using (134) we get i.e.
Remark 7 Returning to (29), we see that namely that Since, by definition, (140) can be written in the comfortable form showing that f (σ ) can be obtained directly from f k (σ ) as limit of the series (142).

A modified CoB for global models
A gravity model, as recalled in the Introduction, is just a finite degree potential (see 1). The reason to concentrate on this particular case is twofold: on the one hand, this is what we really want to get in our analysis, namely that CoB is converging for a global model. On the other hand, we will introduce in the space of global models of degree N a "natural" topology, which however would become more problematic to be extended to the whole L 2 (σ ). So we first define the space of global models of degree N with Then we observe that, given any smooth star-shaped surface, S ≡ {r = R σ , R σ > 0}, u ∈ M N is totally identified by its trace on S, i.e. we have the one-to-one correspondences A rapid inspection of the CoB scheme in Fig. 1 shows that (145) can be formally written as each operator in (146) is invertible, whence the above correspondence.
We call and we notice that M N , H N , H N are all (N +1) 2 dimensional linear spaces, since their elements depend ultimately on the (N + 1) 2 real constants u nm . The one-to-one relations (146) can therefore be used to introduce a norm in one of them and then extend it to the others by isometry. In particular, we will define Such norms correspond to the usual scalar product in L 2 (σ ) such that is an orthonormal basis of H N , and to the unusual scalar product It is clear that in this case is an orthonormal basis of H N . It is clear from (149) and (152) that if we tried to go to infinite dimensions, N → ∞, with this definitions we would end up with functions v defined on S, that in any way can be harmonically continued down to S R ; at this point it is more handy to work directly with finite degree potentials, what justifies our previous choice. Now let us return to the CoB scheme of Fig. 1. The iteration is accomplished by going from S to S R with the pull back operator PB and then going back to S by the upward operator U . Spectrally this means (see 38) As we see from (155) we always have however in general we have that, unless c n q n+1 σ is constant, i.e. S is a sphere, PB sends functions in H N to a subspace of L 2 (σ ), spanned by many more harmonics. The ordinary practice at this point would be that we restrict R σ to be a surface with finite spectrum, i.e. R σ ∈ H p for some p. Yet this implies that q σ = R/R σ has full spectrum up to infinity and the same would be true for f = s u (Bucha et al. 2019). So to cut short the discussion, we make here the unconventional, but numerically equivalent, hypothesis that directly q σ ∈ H p for some p. Then, returning to our discussion, even assuming that we have namely PB : and This means that at every cycle in the iteration we are accumulating more and more high degrees in the residuals which are not useful in building the approximation to f , which we know a priori to have zero components on all the degrees beyond N . An easy remedy to this anomalous inflation of the norm of the residuals is just to truncate U PB f at the maximum degree N . This can be done in two ways, either by projecting U PB f on H N , using the orthonormal basis {q n+1 σ Y nm } and the scalar product ·, · H N , or by projecting PB f on H N , using the orthonormal basis {Y nm } and the scalar product ·, · H N = ·, · L 2 (σ ) . The two results coincide, so we can choose the second approach and introduce the truncation operator T N with kernel to be used with the ordinary L 2 (σ ) scalar product.
We find then, namely the matrix K jk,nm performing the iteration f → K f ∈ H N , is given by Since the scalar product ·, · is coupled with the norm of f in H N , we also have that where the matrix norm is given by the square root of the maximum eigenvalue, i.e. the spectral radius of the matrix.

Remark 8
It is convenient to elaborate already here the form of the matrix K jk,nm when c n is not just 1, but it is given by c n = 1 + (n + 1) h σ R σ . Recalling that h σ = R σ − R and noticing that we easily find c n q n+1 that is the formula we will use in the next section.

Remark 9
It is important to understand that when we pass to finite dimensional potentials and on the same time we restrict the CoB algorithm to a finite dimensional subspace of L 2 (σ ), i.e. H N , we represent the operator K by a finite matrix which has only a finite number of singular values. So the argument that necessarily they tend to 1 used to claim that K = 1 (or larger) (see 100) is not valid here and in fact, if the sandwich conjecture (see 101 and 102) holds, we can expect Indeed, in lack of a rigorous proof, the contractive character of the finite dimensional K expressed in (169) is just an expectation, even if the sandwich conjecture should hold true.

Numerical experiments
In this section, we numerically examine the conjecture (169), the validity of which is critical for CoB to converge. We start by studying the pull back operator (35), for which c n q n+1 σ = q n+1 σ (see 40). The experiments are performed for surfaces S of various complexity: a synthetic topography (Sect. 6.1) and the topographies of the asteroids Bennu (Sect. 6.2) and Eros (Sect. 6.3). Finally, in Sect. 6.4, we briefly report the results obtained with the improved pull back operator (36) with c n q n+1 σ being given by (168). Before moving to the experiments, let us emphasize that numerical investigations of this kind can never prove or disprove the convergence of CoB. Nevertheless, a rigorous analytical study with realistic surfaces S appears to be too difficult, so the numerics is used to shed some light on the behaviour of CoB in practice. Hopefully, the experiments will allow us to draw some general conclusions.

Synthetic topography
Let S be represented by a zonal function as The shape of S is dominated by the oscillations of the Legendre polynomial P M (cos ϑ), while the parameters a and α are numerical constants to control the oscillations. In any case, a and α are chosen such that q σ (ϑ) < 1 for all ϑ (see 11). We call this surface a synthetic topography of degree M.
As an advantage, the synthetic topography is independent of the longitude, allowing us to think of CoB as being independent of the longitude, too. This reduces the size of the matrix K jk,nm and simplifies the evaluation of its elements to (see 165) Note that the (2π) −1 factor is missing before the integrals due to the independence on the longitude. For n = j and q σ = 1, we therefore have 2n + 1 2 π 0 P n (cos ϑ) P n (cos ϑ) sin ϑ dϑ = 1.
The matrix K jk was computed for two synthetic topographies of degree M = 180, one with a 1 = 0.001 and α 1 = 0.0003 and the other with a 2 = 0.01 and α 2 = 0.003. On the Earth, the oscillations of the two topographies would roughly correspond to the range of max σ R σ − min σ R σ ≈ 9 km and 90 km, respectively, with the half-wavelength of 110 km. The larger is the range max σ R σ − min σ R σ , the more non-spherical is the surface S which, in turn, increases the numerical complexity. The {K jn } elements were obtained by (171) for j = 0, 1, . . . , N and n = 0, 1, . . . , N , where N gradually varies from 0 to 720. The scalar products were computed in Fortran in double and quadruple precision using the numerical integration technique of Fukushima (2017). The relative error tolerance δ, affecting the accuracy of the numerical integration, was set to δ = 10 −16 in double precision and to δ = 10 −32 in quadruple precision. In double precision, the norms K were computed for all N = 0, 1, . . . , 720 in Python using the norm function (the ord=2 flag) from the numpy.linalg module (Harris et al. 2020). In quadruple precision, the norms were computed in MATLAB using the ADVANPIX toolbox (www.advanpix.com). Figures 3 and 4 show that somewhere around degree 700, the norms in double precision do not satisfy the (169) conjecture, as they are equal to or larger than the threshold of 1. In quadruple precision, the inequality (169) is, however, satisfied for all studied N . Therefore, we attribute the jump near degree 700 to numerical errors caused by the limited accuracy of the double precision environment. A comparison of the two figures also shows that the norms of the second synthetic topography approach the threshold of 1 faster, owing to the stronger oscillations of S.
As a conclusion, we expect the conjecture (169) is satisfied for the studied synthetic topographies, even for the strongly Fig. 3 Norms K as a function of the maximum harmonic degree of the scalar products N for a synthetic topography (170) defined by M = 180, a 1 = 0.001 and α 1 = 0.0003. The abbreviations dp and qp stand for the double and quadruple precision environments, respectively, in which the scalar products and the norms were computed Fig. 4 The same as Fig. 3 but for a more oscillatory synthetic topography defined by M = 180, a 2 = 0.01 and α 2 = 0.003 oscillating one. This conclusion is valid at least up to the studied maximum degree N = 720.

Asteroid Bennu
In this section, the surface S is represented by a realistic shape model of the asteroid Bennu. We define R σ as where r nm are spherical harmonic coefficients of Bennu's shape taken from the Bucha and Sansò (2021) study. This scenario is significantly more challenging than that with the synthetic topographies for two reasons. First, R σ now varies both with the colatitude and longitude, so that all {K jk,nm } elements up to some maximum degree N are considered. Second, the ratio between the radii of the minimum Brillouin sphere and the maximum Bjerhammar sphere, the shape ratio R/R 0 (see Bucha and Sansò 2021), is about 282.5 m/224 m ≈ 1.261 for Bennu while it is only about (6378 km + 90 km)/6378 km ≈ 1.014 for the second synthetic topography from Sect. 6.1. The shape ratio reflects the extent, to which the surface S deviates from a sphere. Large shape ratio therefore implies highly non-spherical surface S, hence increased numerical complexity (see Sect. 6.1). Due to the complexity of this scenario, we were able to compute K jk,nm only up to degree N = 65. To get the scalar products, we used δ = 10 −8 in the numerical integration technique at first and executed the computation in double precision. Despite the seemingly low maximum degree of 65, the norms were revealed to exceed the threshold of 1 already at degree N = 31 (Fig. 5). To make sure this result is not affected by numerical inaccuracies, we tightened the relative error tolerance to δ = 10 −16 . With the improved numerical accuracy, it turned out the norms are in fact smaller than 1 for all N ≤ 63, meaning that the critical degree moved to higher figures (see Fig. 5). It would certainly be of interest to repeat the experiment in quadruple precision with some δ < 10 −16 to confirm/disprove that the degree of 63 is not again limited by the numerics, similarly as we did in Sect. 6.1. Unfortunately, this is not possible for us with our hardware and software, as already the integration in double precision with δ = 10 −16 required a few weeks of wall-clock time when using about 96 CPU cores. This clearly demonstrates that our numerical experiments cannot provide definite answers on the validity of (169), but at most some hints, though valuable.
To summarize the experiment, we have shown that even with as irregularly shaped body as the asteroid Bennu, the norms stay below the threshold of 1 at least up to degree 63. Using the very same shape of Bennu, Bucha and Sansò (2021) have shown that a relative accuracy of ∼10 −6 can be achieved on S in terms of the gravitational potential that is expanded in a single series of external spherical harmonics up to degree 65. This might indicate that CoB could deliver reasonably accurate gravity field models even in advanced conditions.

Asteroid Eros
Finally, we repeated the previous experiment, but this time with a non-spherical asteroid Eros. The shape coefficients r nm of Eros were taken from the Zuber et al. (2000) study up to degree 5. The shape ratio for such defined surface reads roughly 16000 m/3250 m ≈ 4.923 which demonstrates the non-sphericity of Eros.
The K jk,nm elements were computed up to degree N = 20 in double precision with δ = 10 −16 and in quadruple precision with δ = 10 −32 . In both cases, the norms are almost identical and quickly approach the threshold of 1, exceeding it already at degree 6 ( Fig. 6). This indicates that there might indeed exist conditions, for which the conjecture (169) is not satisfied. Yet, we cannot reject the conclusion that (169) is satisfied but instead quadruple precision and the relative error tolerance of δ = 10 −32 still do not ensure sufficient numerical accuracy for Eros.

Improved pull back operator
To study the effect of the Taylor order on the downward continuation, we repeated some of the previous experiments with the improved pull back operator (36). We studied the norms for the second synthetic topography in quadruple precision (see Sect. 6.1) and for Eros in double precision (Sect. 6.3). The plots of the norms are not shown here as the curves are visually the same as seen in the previous sections.
With the second synthetic topography, the norms stayed below the threshold of 1 with both variants of the pull back operator. However, the improved pull back operator produced slightly smaller norms (e.g., 0.999 vs. 0.994 at degree 720 in favour of the improved pull back operator). In case of Eros, the norms exceeded the threshold of 1 at degree 7 which is one degree higher than with the pull back operator (35) (see Sect. 6.3). The generally smaller norms indicate that the improved pull back operator (36) has better convergence properties than (35), though the enhancement is probably not significant.

Some conclusions and perspectives
In the paper the concept of change of boundary has been revisited and in particular the functional implications of the hypothesis that have been explored, showing that even in the case of a nonexpansive operator, K = 1, we can have a convergence of the method when the concept of Cesaro limit is used. The hypothesis (174) has no strict analytical proof for the moment, yet it is a consequence of a first-order perturbative argument, in which the surface S ≡ {r = R σ } is supposed to be not too far from a sphere.
In lack of a rigorous proof, our hypothesis has been subject to numerical tests in order to verify its plausibility.
In the first experiment, we have used data from a simulation of the gravity field of objects with various topographic surfaces with a width up to 1.4 % of the radius (∼90 km in the case of the Earth) and wavelength corresponding to a maximum degree 720.
In the second experiment, data related to the asteroid Bennu (Bucha and Sansò 2021) have been used. In this case the topographic surface is quite smooth but the ratio between minimum Brillouin radius and maximum Bjerhammar radius, R/R 0 , is about 1.261 meaning a significant, but not too large, deviation from a spherical shape. In both experiments, the hypothesis (174) has proved to be tenable.
In a third experiment we have been stressing the conditions of validity of hypothesis (174) by using data related to the asteroid Eros which has a shape very far from a sphere, with ratio R/R 0 ≈ 4.923. In this case the K norm has been estimated to be larger than 1, i.e. the CoB theory is not applicable.
Furthermore two forms of the pull back operator have been tested, arriving at the conclusion that they display a not significantly different behaviour as for what concerns convergence.
Summarizing we may claim that the hypothesis (174), and the consequent convergence of the method, can be considered as effective when working with boundaries S not far from the sphere, as for long wavelength behaviours, while short wavelength perturbations seem not to effect too much convergence.
On the contrary, for a surface S significantly deformed with respect to the sphere, we certainly cannot rely on the CoB method and we expect non-convergence even for lowdegree truncated models.
The threshold between the two behaviours still needs to be studied both theoretically and numerically. Acknowledgements The computations were performed at the HPC centre at the Slovak University of Technology in Bratislava, which is a part of the Slovak Infrastructure of High Performance Computing (SIVVP Project, ITMS code 26230120002, funded by the European region development funds, ERDF). The plots were prepared with the Python's Matplotlib module (Hunter 2007).

Author Contributions
Although FS has worked more on the theoretical development of the study and BB on the numerical experiments, a continuous discussion between the authors has made the paper an equal result of theirs.
Funding Open access funding provided by The Ministry of Education, Science, Research and Sport of the Slovak Republic in cooperation with Centre for Scientific and Technical Information of the Slovak Republic. BB was supported by the VEGA 1/0809/21 Project.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.