manuscript No. (will be inserted by the editor) Two-parameter regularization of ill-posed spherical pseudo-differential equations in the space of continuous functions Hui Cao · Sergei V. Pereverzyev · Ian H. Sloan · Pavlo Tkachenko 5 Received:date/Accepted:date 1 0 2 Abstract In this paper, a two-step regularization method is used to solve an n ill-posed spherical pseudo-differential equation in the presence of noisy data. a For the first step of regularization we approximate the data by means of a J spherical polynomial that minimizes a functional with a penalty term consist- 2 ing of the squared norm in a Sobolev space. The second step is a regularized ] collocation method. An error bound is obtained in the uniform norm, which A ispotentiallysmallerthanthatforeitherthenoisereductionaloneorthereg- N ularized collocation alone. We discuss an a posteriori parameter choice, and . presentsomenumericalexperiments,whichsupporttheclaimedsuperiorityof h the two-step method. t a m Keywords Two-parameter regularization · spherical pseudo-differential [ equations · quasi-optimality criterion 1 v 2 H.Cao 6 Guangdong Province Key Laboratory of Computational Science, Sun Yat-sen University, 3 Guangzhou,China 0 E-mail:[email protected] 0 . S.V.Pereverzyev 1 Johann Radon Institute for Computational and Applied Mathematics, Austrian Academy 0 ofSciences,Altenbergerstrasse69,4040Linz,Austria 5 E-mail: [email protected] 1 I.H.Sloan : v SchoolofMathematicsandStatistics,UniversityofNewSouthWales,SydneyNSW2052, i Australia X E-mail:[email protected] r P.Tkachenko a Johann Radon Institute for Computational and Applied Mathematics, Austrian Academy ofSciences,Altenbergerstrasse69,4040Linz,Austria E-mail:[email protected] 2 HuiCaoetal. 1 Introduction Mathematical models appearing in geoscience commonly have the form of an ill-posed spherical pseudo-differential equation Ax=y, (1) where A is a pseudo-differential operator that relates continuous functions x∈C(Ω ) and y ∈C(Ω ) defined on concentric spheres Ω ,Ω ∈R3 of radii R ρ R ρ R ≤ ρ. For example, in satellite geodesy, this approach has been introduced in [5,25] where the spheres Ω ,Ω are models of respectively the surface of R ρ the Earth and the surface traversed by the satellite. Recall that a spherical pseudo-differential operator A:C(Ω )→C(Ω ) is R ρ a linear operator that assigns to any x∈C(Ω ) a function R ∞ 2k+1 (cid:18) (cid:19) X X 1 · Ax:= a x Y ∈C(Ω ), (2) k bk,jρ k,j ρ ρ k=0 j=1 where (cid:28)1 (cid:16) · (cid:17) (cid:29) 1 Z (cid:16)τ (cid:17) x = Y ,x(·) := x(τ)Y dΩ (τ) bk,j R k,j R R k,j R R L2(ΩR) ΩR are the spherical Fourier coefficients, and Y (·),j = 1,2,...,2k+1, are the k,j sphericalharmonics[16]ofdegreek whichareL -orthonormalwithrespectto 2 the unit sphere Ω ∈R3, as a result of which 1 (cid:28)1 (cid:16) · (cid:17) 1 (cid:16) · (cid:17)(cid:29) Y , Y =hY ,Y i =δ δ . (3) R k,j R R k0,j0 R k,j k0,j0 L2(Ω1) kk0 jj0 L2(ΩR) Thesequenceofrealnumbers(a )∞ isreferredtoasthesphericalsymbolof k k=0 A. We shall assume that a is positive, and converges monotonically to zero. k In the case when the symbol sequence a tends to zero fast enough the k operator A is compact. Therefore, its inverse A−1 is unbounded, and keeping inmindHadamard’sdefinitionofawell-posedproblem(existence,uniqueness, andcontinuityofinverse),weconcludethatforthiscasethefirstandthethird conditions are violated, and the problem (1), (2) with y ∈C(Ω ) becomes ill- R posed. Therefore, a regularization technique should be employed for solving it [4]. As examples of the ill-posed problem (1), (2) we can mention the satellite- (cid:16) (cid:17)k to-satellitetrackingproblem(SST-problem)witha = k+1 R ,thesatellite k ρ ρ (cid:16) (cid:17)k gravity gradiometry problem (SGG-problem) with a = (k+1)(k+2) R , etc. k ρ2 ρ (formoredetailsontheseandotherexampleswecanreferaninterestedreader to [5]). These problems are severely ill-posed because of the occurrence of the geometric factor (R/ρ)k. Two-parameterregularization 3 Itisworthmentioningthatmanypracticalapplicationsuseafinitedimen- sional approximation of the solution of (1), (2). For example, Earth Gravity Models, such as EGM96 or EGM2008 [20] are parametrized by the spherical Fourier coefficients up to some prescribed degree M. Note that in applications the function y is assumed to be continuous. However, in practice one is provided just with a finite number of points {t }N ⊂Ω atwhichinformationaboutthevaluesofy iscollected.Itshould i i=1 ρ benotedalsothatthepointwisedatay(cid:15)(t )containmeasurementerrors,which i can be modeled, for example, in the following way: |y(cid:15)(t )−y(t )|≤(cid:15) , i=1,...,N. i i i We assume for convenience that there exists a function y(cid:15) ∈ C(Ω ) standing ρ for the noisy version of the original function y, such that ky(cid:15)−yk ≤(cid:15):= max {(cid:15) }, C(Ωρ) 1≤i≤N i where (cid:15) are measurement errors. i In such a setup the problem (1), (2) is reduced to the following spherical pseudo-differential operator equation M 2k+1 (cid:18) (cid:19) X X 1 · A x:= a x Y =y(cid:15). (4) M k bk,jρ k,j ρ k=0 j=1 Severalregularizationtechniquescanbeusedfortreating(4).Forexample, the method where the discretization level M plays the role of regularization parameter was discussed in [2,3,9]. However, in our case we assume that the value of M is prescribed, thus this approach cannot be used. The problem (4), but without noise in the right-hand side, was studied in [13]. In this noise-free case the solution x can be approximated by solving the equation AMx=VbMy, where VbM : C(Ωρ) → PM(Ωρ) is the so-called quasi-interpolatory operator. Here by P (Ω ) we denote the set of all spherical polynomials of degree less M ρ thanorequaltoM,orinotherwordstherestrictiontoΩ ofthepolynomials ρ in R3 of degree less than or equal to M. Thus, in [13] the authors suggest constructingapolynomialapproximationofy fromtheoriginalpointwisedata y(t ),i=1,2,...,N, and then formally inverting the operator A . i M In principle, this idea can be used also for the ill-posed case (4). However, in contrast to the approximation of a noise-free continuous function y on the sphere by means of polynomials, which has been discussed by many authors (see,forexample,[7,24,26,28,30]),theapproximationofnoisyfunctionsy(cid:15) has been studied only recently in [1,21], to the best of our knowledge. Applying the method from [1,21] we first perform so-called data-smoothing (or noise reduction). After this data preprocessing the formal inversion of A should M 4 HuiCaoetal. be safer. However, for performing the noise reduction step one needs a priori information about the smoothness of the function y, which is usually not available,afactthatmakesthedirectapplicationofthescheme[21]notalways appropriate. In a different direction, for estimating the Fourier coefficients x directly bk,j from noisy measurements y(cid:15)(t ) a regularized collocation method has been re- i centlypresentedin[18].Thismethodisbasedonthestandardandwidelyused Tikhonov-Phillipsregularization.However,itiswell-knownthattheTikhonov- Phillips method suffers from saturation [15,19], meaning that the accuracy of reconstruction cannot be improved regardless of the smoothness of the solu- tion x. Another point is that while the Tikhonov-Phillips method has been well studied in the Hilbert space L , to the best of our knowledge, no analysis 2 has been done in the space of continuous functions, a natural choice for our noise model. In the present study we combine these two approaches. Moreover, in con- trasttothepreviousresults,wewillanalyzetheapproximationinthespaceof continuous functions. A combination of two regularization methods into one can be seen as a two-step regularization, in which we use the composition R ◦T of data smoothing operator T , and a regularized collocation α,M λ,M λ,M operator R . In the literature there are not so many studies on two-step α,M regularization,andwecanreferonlyto[11,12].However,theanalysisinthose papers does not correspond to the setting of our problem (1)–(4). The paper is organized as follows. In the next section we present the regu- larization method for noise reduction and define the data smoothing operator T . In Section 3 we will give a short overview of the regularized collocation λ,M method,anddefinetheoperatorR .Section4isdevotedtotheoreticalerror α,M bounds for the constructed two-parameter regularization. We will show that the approximation has the potential to perform at least as well as the better of the one-parameter regularizations which are involved in the composition. Finally, in the last section we discuss an a posteriori parameter choice rule, and present some numerical experiments supporting the claimed superiority of our method. 2 Data noise reduction At the first step of our scheme we approximate the noisy continuous function y(cid:15) ∈C(Ω )bymeansofasphericalpolynomialp ∈P (Ω ).Asdiscussedin ρ M M ρ the Introduction, instead of y we are provided only with pointwise measure- ments y(cid:15)(t ),i = 1,2,...,N. Therefore we introduce the sampling operator i S :C(Ω )→RN for which N ρ ω S y(cid:15) :=(y(cid:15)(t ),y(cid:15)(t ),...,y(cid:15)(t )). N 1 2 N By RN we denote the vector space RN equipped with the inner product ω Two-parameterregularization 5 N X hη,γi := ω η γ , η,γ ∈RN; ω i i i i=1 the corresponding norm is kηk = hη,ηi1/2. Here ω ,ω ,...,ω are positive ω ω 1 2 N weights in a cubature formula, and t ,t ,...,t ∈ Ω are the corresponding 1 2 N ρ cubaturepoints,withthecubaturerulehavingthepropertyofbeingexactfor all polynomials of degree up to 2M, N Z X ∀p∈P (Ω ), ω p(t )= p(ζ)dΩ (ζ). (5) 2M ρ i i ρ i=1 Ωρ A method for approximating y from S y(cid:15) ∈ RN based on such a cubature N ω rule (but with all weights equal) was proposed recently in [1], while its gener- alization to other positive weights was presented in [21]. We briefly outline the method from [21] and its extension to pseudo- differential operator equations. We start with the observation that the space P (Ω ) of spherical polynomials p of degree at most M can be considered as M ρ the reproducing kernel Hilbert space (RKHS) H generated by the kernel K M 2k+1 (cid:18) (cid:19) (cid:18) (cid:19) X X 1 t 1 τ K(t,τ)= β−2 Y Y , t,τ ∈Ω , (6) k ρ k,j ρ ρ k,j ρ ρ k=0 j=1 where β = (β ,β ,...,β ,...) is a non-decreasing sequence of positive pa- 1 2 M rameters. Note that the inner product of H associated with this kernel can K be written as XM β2 2Xk+1(cid:28) (cid:18)·(cid:19) (cid:29) (cid:28) (cid:18)·(cid:19) (cid:29) hf,gi = k Y ,f Y ,g . (7) HK ρ2 k,j ρ k,j ρ k=0 j=1 L2(Ωρ) L2(Ωρ) The reproducing property hf,K(·,τ)i = f(τ) for τ ∈ Ω can be verified HK R easily. In this paper the numbers β ,β ,...,β are assumed to be given a 1 2 N priori. We shall see their role later. The approximation p studied in [1,21] appears as the minimizer of the min following functional n o p =argmin kS p−S y(cid:15)k2 +λkpk2 , p∈P (Ω ) , (8) min N N ω HK M ρ where λ ≥ 0 is a regularization parameter. Note that (8) can be seen as the definition of the data smoothing operator T :C(Ω )→P (Ω ) such that λ,M ρ M ρ p =T y(cid:15). Note that the regularization term involves the H -norm of p, min λ,M K and hence depends on the rate of growth of the parameters β . The solution k of (8) is given by the following theorem. 6 HuiCaoetal. Theorem 1 Assume that the points {t } and weights {ω } are such that (5) i i holds. Then the minimizer T y(cid:15) =p in (8) has the form λ,M min M 2k+1 (cid:18) (cid:19) N (cid:18) (cid:19) (T y(cid:15))(·)=X 1 X 1Y · Xω 1Y ti y(cid:15)(t ). (9) λ,M 1+λβ2 ρ k,j ρ iρ k,j ρ i k=0 k j=1 i=1 Proof It is known that the minimizer of the functional in (8) can be written as T y(cid:15) =(λI+S∗S )−1S∗S y(cid:15), (10) λ,M N N N N where S∗ : RN → H is the adjoint of S , and I is the identity operator in N ω K N H . By definition, for (η )N ∈RN we have K i i=1 ω N N X X hη,S fi = η f(t )ω = η hK(t ,·),f(·)i ω N ω i i i i i HK i i=1 i=1 * N + X = ω K(t ,·)η ,f(·) , i i i i=1 HK thus there holds N X (S∗η)(·)= ω K(t ,·)η , ∀η ∈RN, N i i i ω i=1 and thereby N X (S∗S f)(·)= ω K(t ,·)f(t ). N N i i i i=1 Inserting this expression into (10), we observe that p solves min N N X X λp (·)+ ω K(t ,·)p (t )= ω K(t ,·)y(cid:15)(t ) (11) min i i min i i i i i=1 i=1 The spherical harmonic expansion of the polynomial p =T y(cid:15) is min λ,M M 2k+1(cid:28) (cid:18) (cid:19)(cid:29) (cid:18) (cid:19) X X 1 · 1 t p (t)= p , Y Y , t∈Ω . (12) min min ρ k,j ρ ρ k,j ρ ρ k=0 j=1 L2(Ωρ) To find the coefficients in this expansion, we write the second term on the left-hand side of (11), on using (5) and then (6), as N X ω K(t ,t)p (t )=hp ,K(·,t)i i i min i min L2(Ωρ) i=1 N 2k+1(cid:28) (cid:18) (cid:19)(cid:29) (cid:18) (cid:19) X X 1 · 1 t = β−2 p , Y Y . k min ρ k,j ρ ρ k,j ρ k=0 j=1 L2(Ωρ) Two-parameterregularization 7 After expanding the right-hand side of (11) using (6), and then equating the (cid:16) (cid:17) coefficients of Y t , we obtain k,j ρ (cid:28) (cid:18) (cid:19)(cid:29) N (cid:18) (cid:19) t X t (λ+β−2) p ,Y =β−2 ω Y y(cid:15)(t ), t∈Ω , k min k,j ρ k i k,j ρ i ρ L2(Ωρ) i=1 which on substituting into (12) yields the desired result. tu At this point the solution of the first step depends on the regularization parameter λ, and the penalization weights β . As mentioned in the Introduc- k tion, the choice of the regularization parameter λ will be addressed in the last section. A data-driven choice of the penalization weights β has been recently k discussed in [21]. The assumption that the function x on Ω is continuous implies that R x∈L2(Ω ),andhencethatitsFouriercoefficients(cid:10)1Y (cid:0) · (cid:1),x(cid:11) with R R k,j R L2(ΩR) respect to the basis of spherical harmonics are square-summable, i.e, X∞ 2Xk+1(cid:12)(cid:12)(cid:28)1 (cid:16) · (cid:17) (cid:29) (cid:12)(cid:12)2 (cid:12) Y ,x (cid:12) <∞. (cid:12) R k,j R (cid:12) k=0 j=1 (cid:12) L2(ΩR)(cid:12) Any additional smoothness of x can be measured in terms of the summability ofFouriercoefficientswithsomeincreasingweightsdependingonthesequence (β ) or, as it is usual for the regularization theory (see, e. g., [14]), on the k symbol(a ).Therefore,itisconvenienttointroducetwoSobolevspacesWφ,β, k and Wψ,a, the first depending on the sequence (β ) and the second on the k symbol (a ), and defined by k Wφ,β :=g∈L2(ΩR):Xk∞=02Xjk=+11(cid:12)(cid:12)(cid:12)(cid:10)R1Yk,j(cid:0)φR2·((cid:1)βk,−g2(cid:11))L2(ΩR)(cid:12)(cid:12)(cid:12)2 =:kgk2Wφ,β <∞, (13) Wψ,a:=g∈L2(ΩR):X∞ 2Xk+1(cid:12)(cid:12)(cid:12)(cid:10)R1Yk,j(cid:0)ψR·2(cid:1)(a,2g)(cid:11)L2(ΩR)(cid:12)(cid:12)(cid:12)2 =:kgk2Wψ,a <∞,(14) k=0 j=1 k where φ,ψ are non-decreasing functions such that φ(0) = 0 and ψ(0) = 0. In the literature, see, e.g., [14], the functions φ,ψ go under the name of “smoothness index functions”, and are usually unknown. 3 Regularized collocation method after noise reduction After the first step of our method the reduced original equation (4) is trans- formed into the equation 8 HuiCaoetal. A x=T y(cid:15). (15) M λ,M The regularized collocation method (see [18]) is now applied to this equation, yielding an approximate solution x(cid:15) of the equation (15), defined as the α,λ minimizer of the functional kB x−S T y(cid:15)k2 +αkxk2 , (16) N,M N λ,M ω L2(ΩR) where B := S A : L (Ω ) → RN, and α ≥ 0 is the regularization N,M N M 2 R ω parameter for the second step. The minimizer of (16) can be written in the form x(cid:15) =(αI+B∗ B )−1B∗ S T y(cid:15), (17) α,λ N,M N,M N,M N λ,M where B∗ :RN →L (Ω ) is the adjoint of B , given by N,M ω 2 R N,M (B∗ η)(·)=XM a 2Xk+1 1Y (cid:16) · (cid:17)XN ω 1Y (cid:18)ti(cid:19)η , N,M k R k,j R iρ k,j ρ i k=0 j=1 i=1 from which it follows, using (4), that M 2k+1 X X 1 (cid:16) · (cid:17) B∗ B x= a2 Y x . N,M N,M k R k,j R bk,j k=0 j=1 Using (17) and (9) we then obtain explicitly x(cid:15) :=XM ak 1 2Xk+1 1Y (cid:16) · (cid:17)XN ω 1Y (cid:18)ti(cid:19)y(cid:15)(t ). (18) α,λ α+a2 1+λβ2 R k,j R iρ k,j ρ i k=0 k k j=1 i=1 If we now define R :P (Ω )→P (Ω ) as α,M M ρ M R (R p)(·):=(cid:16)(cid:0)αI+B∗ B (cid:1)−1B∗ S p(cid:17)(·) α,M N,M N,M N,M N = XM ak 2Xk+1 1Y (cid:16) · (cid:17)XN ω 1Y (cid:18)ti(cid:19)p(t ) (19) α+a2 R k,j R iρ k,j ρ i k=0 k j=1 i=1 = XM ak 2Xk+1 1Y (cid:16) · (cid:17)(cid:28)1Y ,p(cid:29) . (20) α+a2 R k,j R ρ k,j k=0 k j=1 L2(Ωρ) then we can write x(cid:15) = R T y(cid:15). Thus R reflects the regularized α,λ α,M λ,M α,M collocation second step of the method, whereas T corresponds to the noise λ,M reduction first step. Asonecanseefrom(18),intheendweapproximatetheunknownsolution x on Ω by means of a spherical polynomial of degree M. In this context a R question arises about the best approximation of a continuous function x by means of a polynomial of degree M. In turn the quality of best polynomial Two-parameterregularization 9 approximation is determined by the smoothness of x. We shall assume that x liesintheintersectionofWφ,β andWψ,a,see(13)or(14),thusthesmoothness ofxisencodedinφandβ ,orψ anda .Forexample,ifthesmoothnessindex k k function φ(t) and the sequence β ={β } increase polynomially with t and k, k then Jackson’s theorem on the sphere (see [22], Theorem 3.3) tells us that for x∈Wφ,β there is ν >0 such that inf kx−pk =O(cid:0)M−ν(cid:1). (21) p∈PM C(ΩR) At the same time, if the sequence β = {β } increases exponentially then k for polynomially increasing φ and x∈Wφ,β we have inf kx−pk =O(cid:0)e−qM(cid:1), p∈PM C(ΩR) where q is some positive number that does not depend on M. Intheerroranalysisofthenextsectionweshallmakeuseofaconstructive polynomial approximation introduced in [26], in which x∈C(Ω ) is approxi- R mated by V x∈P (Ω ), given by M M R XM (cid:18) k (cid:19)2Xk+1 1 (cid:18) t (cid:19)(cid:28)1 (cid:16) · (cid:17) (cid:29) (V x)(t)= h Y Y ,x , t∈Ω , M M R k,j R R k,j R R k=0 j=1 L2(ΩR) (22) where h (a “filter function”) is a continuously differentiable function on R+ satisfying (cid:26) 1, t∈[0,1/2], h(t)= 0≤h(t)≤1fort∈R+ 0, t∈[1,∞), Explicit examples of suitable filter functions h can be found in [27]. It is im- portantforourlateranalysisthat,asshownin[26],V xis(uptoaconstant) M an optimal approximation in the uniform norm, in the sense that, kx−V xk ≤c inf kx−pk , (23) M C(ΩR) p∈P[M/2] C(ΩR) where c is a generic constant, which may take different values at different occurrences, and [·] denotes the floor function. In view of (21), for polynomially increasing φ,β and for x∈Wφ,β we have kx−V xk ≤c[M/2]−ν ≤cM−ν. M C(ΩR) Ontheotherhand,forexponentiallyincreasingβ andpolynomiallyincreasing φthetheory[23]suggeststakingh(t)≡1,t∈[0,1].Inthiscasetheright-hand √ sideof (23)hastobemodifiedbymultiplyingby M,andbyreplacing[M/2] by M, thus for x∈Wφ,β there holds √ √ kx−V xk ≤c M inf kx−pk ≤c Me−qM. M C(ΩR) p∈PM C(ΩR) 10 HuiCaoetal. 4 Error estimation in uniform norm Inthissectionwewillestimatetheuniformerrorofapproximationofxbythe polynomial x(cid:15) given by (18). It is clear that α,λ (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)x−x(cid:15)α,λ(cid:13)C(ΩR) =(cid:13)x−Rα,MTλ,My(cid:15)(cid:13)C(ΩR) (24) (cid:13) (cid:13) ≤kx−VMxkC(ΩR)+(cid:13)VMx−Rα,MUMy(cid:13)C(ΩR) (cid:13) (cid:13) +(cid:13)(Rα,M −Rα,MTλ,M)UMy(cid:13) C(ΩR) (cid:13) (cid:13) (cid:0) (cid:1) +(cid:13)Rα,MTλ,M(cid:13)C(Ωρ)→C(ΩR) kUMy−ykC(Ωρ)+ky−y(cid:15)kC(Ωρ) , where U y is the same as V x in (22) except that it is an approximation on M M the larger sphere Ω instead of on Ω , ρ R M (cid:18) (cid:19)2k+1 (cid:18) (cid:19)(cid:28) (cid:18) (cid:19) (cid:29) X k X 1 t · (U y)(t)= h Y Y ,y(·) . (25) M M ρ2 k,j ρ k,j ρ k=0 j=1 L2(Ωρ) It is natural to assume that kU y−yk < (cid:15), since otherwise data noise M C(Ωρ) is dominated by the approximation error and no regularization is required. We also restrict ourselves to the case when kx−V xk < c(cid:15), otherwise M C(ΩR) the term kx−V xk , representing the error of almost best polynomial M C(ΩR) approximation, will dominate the error bound. Then the bound (24) can be reduced to the following one (cid:13)(cid:13)x−x(cid:15)α,λ(cid:13)(cid:13)C(ΩR) ≤kVMx−Rα,MUMykC(ΩR) (26) +k(R −R T )U yk +c(1+kR T k )(cid:15). α,M α,M λ,M M C(ΩR) α,M λ,M C(Ωρ)→C(ΩR) An estimate of the coefficient in the last term is given by the following theorem. Theorem 2 Under the conditions of Theorem 1 (cid:12) (cid:12) kRα,MTλ,MkC(Ωρ)→C(ΩR)≤ R1ρtm∈ΩaxR(cid:12)(cid:12)(cid:12)(cid:12)XN ωiXM 4π(α(+2ka+2k)(11)a+kλβk2)Pk(cid:16)tR·ρti(cid:17)(cid:12)(cid:12)(cid:12)(cid:12), i=1 k=0 where P are the Legendre polynomials of degree k. k Proof Inviewofthedefinition(9)and(19),togetherwith(3)withRreplaced by ρ and the property (5) of the cubature rule, we can write (R T y)(t) α,M λ,M M 2k+1 (cid:18) (cid:19) N (cid:18) (cid:19) = X ak 1 X 1Y t Xω 1Y ti y(t ) α+a2 1+λβ2 R k,j R iρ k,j ρ i k=0 k k j=1 i=1 M N (cid:18) (cid:19) = 1 X2k+1 ak 1 Xω P t·ti y(t ), ∀t∈Ω , Rρ 4π α+a2 1+λβ2 i k Rρ i R k=0 k k i=1