Approximate solution of two dimensional disc-like systems by one dimensional reduction: an approach through the Green function formalism using the Finite Elements Method
aa r X i v : . [ phy s i c s . c o m p - ph ] S e p Approximate solution of two dimensional disc-like systems by one dimensionalreduction: an approach through the Green function formalism using the FiniteElements Method
Alejandro Ferrero Botero ∗ Departamento de Ciencias B´asicas, Universidad Cat´olica de Colombia - Bogot´a, Colombia
Juan Pablo Mallarino † Facultad de Ciencias — Laboratorio Computacional HPC, Universidad de los Andes - Bogot´a, Colombia (Dated: October 6, 2020)We present a comprehensive study for common second order PDE’s in two dimensional disc-like systems and show how their solution can be approximated by finding the Green function ofan effective one dimensional system. After elaborating on the formalism, we propose to securean exact solution via a Fourier expansion of the Green function, which entails to solve an infinitelycountable system of differential equations for the Green-Fourier modes that in the simplest case yieldsthe source-free Green distribution. We present results on non separable systems—or such whosesolution cannot be obtained by the usual variable separation technique—on both annulus and discgeometries, and show how the resulting one dimensional Fourier modes potentially generate a near-exact solution. Numerical solutions will be obtained via finite differentiation using FDM or FEMwith the three-point stencil approximation to derivatives. Comparing to known exact solutions, ourresults achieve an estimated numerical relative error below 10 − . I. INTRODUCTION
In the present work we elaborate on the FEM for solv-ing complex two-dimensional partial differential equa-tions (DE) using a Green function construction. Green’smethod has been employed extensively in Physics forsolving Laplace’s equation and associates in a cornu-copia of areas, such as Quantum and Statistical Mechan-ics. In quantum mechanics, for example, the method ofnonequilibrium Green’s functions (NEGF) has been usedto study the Brownian motion of a quantum oscillator[1], quantum thermal transport [2, 3], derive quantumkinetic equations [4], study hadronic physics [5], amongothers. In statistical mechanics, some of the applicationsof the Green functions include the predictions of some ob-servables [6], help to describe 1D hydrodynamic models[7], finding electrical properties of some physical systems[8, 9], study nonextensive statistical mechanics with newnormalized q -expectation values [10], and so much more.Even the Green functions are used in quantum field the-ory to describe the propagators of quantum fields in theperturbative regime.Not only are Green functions useful to solve systemsdescribed by inhomogeneous differential equations, butthey can also be used to describe thermodynamic prop-erties. For instance, the density and correlations of par-ticles immersed in two dimensional two component plas-mas at certain temperatures can be described by sets ofGreen functions [11–13].In order to study how an inhomogenous partial DE canbe solved by the method we propose, we start defining a ∗ [email protected] † [email protected] differential operatorˆ L { r } (cid:3) = ( ~ ∇ { r } + ~f ( r )) · ( ~ ∇ { r } (cid:3) ) + g ( r ) (cid:3) , (1)acting on a scalar field in ℜ d , with d the dimension of thesystem— i.e. r ∈ ℜ d . This operator is known in othercontexts as the Liouville operator; via this definition, weoften describe the evolution of a relevant quantity ψ bymeans of the equation ∂ t ψ ( r , t ) − ˆ L ψ ( r , t ) = 0 as it isthe case of the wave function in quantum mechanics. Forexample, in the diffusion phenomenon the functions takethe form ~f ( r , t ) = ~ ∇ D ( r , t ) and g ( r ) = 0, and for theHelmholtz equation ~f ( r , t ) = 0 and g ( r ) = m , with m aconstant.Finding solutions to the latter has motivated the devel-opment of numerical methods that grow in number andcomplexity. For instance, using Restricted BoltzmannMachines we can engineer an artificial neural networkthat is able to accurately sample the probability distri-bution for quantum statistical systems [14, 15]. However,some effort can be made from a mathematical point ofview prior to implementing a full scale numerical calcu-lation.Green’s function—or more precisely, distribution—isperhaps the most interesting artifact of a huge bag oftricks that we have when facing differential equations. Itspower relies on the possibility of inverting the differentialoperator ˆ L to solve the inhomogeneous equationˆ L ψ ( r ) = φ ( r ) , (2)with ψ ( r ) and φ ( r ) two scalar functions. Hinting thatits existence, the Green Distribution, is conditioned bysome properties of ˆ L .A disadvantage of the Green methodology is the du-plication of degrees of freedom, encouraging researchersto find ψ ( r ) directly. Our aim is not to develop a gen-eralized theory for an arbitrary problem and number ofdimensions. Despite this, we can look into the conse-quences of breaking down one dimension by focusing onthe simple two-dimensional case.Two dimensional systems are of great interest in statis-tical mechanics [11–13], material sciences [16–18], quan-tum computing [19, 20], high energy physics [21, 22], ionicfluids [23], theoretical mathematics [24, 25], and manyothers.The outline of the paper is as follows. We first remindsome relevant known results for the Green’s function con-struction in section II prior to presenting the strategy tomove from 2D to 1D in section II D. We lay out a clevergeometric interpretation of the result in section II A fol-lowed by a connection to a relevant theory for HilbertSpace functions in section II C. Consequently, we nextdiscuss its implications towards finding the Green func-tion using FDM in section III. We present some mathe-matical results that include the solution of some knownresults for testing purposes, the implementation of themethod in a non-separable 2D system, and a discussionof how the algorithm can be adapted to solve the heatdiffusion problem in thermal equilibrium in section IV.Finally, we wrap up the conclusions in section V. Inter-mediate calculations and numerical details are left forfurther inspection in appendices. II. FRAMEWORK
We start studying the Green function formalism bypostulating the convolution identity from the Dirac dis-tribution, ψ ( r ) = Z r ′ ψ ( r ′ ) δ ( r ′ − r ) w ( r ′ , r ) d r ′ , (3)with w ( r ′ , r ) a weight function properly defined by twoconditions; the first of which w ( r , r ) = 1. Now by defin-ing G ( r ′ , r ) as, ˆ L { r ′ } G ( r ′ , r ) = δ ( r ′ − r ) , (4)with δ ( r ′ − r ) = δ ( r − r ′ ) the ℜ d Dirac delta distribution,then, ψ ( r ) = Z r ′ ψ ( r ′ ) h ˆ L { r ′ } G ( r ′ , r ) i w ( r ′ , r ) d r ′ . (5)The second condition over w ( r , r ′ ) will be determinedin such a way that ˆ L is self-adjoint (Hermitian), or equiv-alently Z r ′ ψ ( r ′ ) h ˆ L { r ′ } G ( r ′ , r ) i w ( r ′ , r ) d r ′ = Z r ′ h ˆ L { r ′ } ψ ( r ′ ) i G ( r ′ , r ) w ( r ′ , r ) d r ′ + b.c. , with added Dirichlet or Neumann boundary conditions(b.c.). Direct substitution into eq. (5), using eq. (2),yields ψ ( r ) = Z r ′ G ( r ′ , r ) φ ( r ′ ) w ( r ′ , r ) d r ′ + b.c. . (6) A. On the nature of w ( r , r ′ ) and ˆ L − This former known result deserves a more delicate look,particularly, on the existence of the weight function, andhow previous solution relates with the usual convolutiontheorem ψ ( r ) = R r ′ G ( r , r ′ ) φ ( r ′ ) d r ′ + b.c.. As mentioned,an appropriate choice for the weight function ensuresthat eq. (5) reproduces eq. (6). This is done by usingthe Green’s and Divergence theorem in eq. (5) to per-form an integration by parts. After simplifications, werealize that by choosing the weight function such that ~ ∇ { r ′ } w ( r ′ , r ) − w ( r ′ , r ) ~f ( r ′ ) = 0 (See section A for fur-ther details) we ensure that the operator is self-adjoint!This is essential to Green’s method. Hence, if no weightfunction exists, we might be forced to use other analyticaland/or numerical procedures in order to find ψ . For thatmatter, the range of problems that we aim to analyze isnarrowed down to the few ones satisfying the aforemen-tioned condition; despite this, a great many of this subsetare of special interest for Mathematics and Physics.Assuming w ( r ′ , r ) exists we are able to incorporate thepremise for eq. (6) yielding exactly, ψ ( r ) = Z r ′ G ( r ′ , r ) φ ( r ′ ) w ( r ′ , r ) d r ′ + I ∂ r ′ w ( r ′ , r ) × h ψ ( r ′ ) ~ ∇ { r ′ } G ( r ′ , r ) − G ( r ′ , r ) ~ ∇ { r ′ } ψ ( r ′ ) i · n d S ′ . (7)Note how the second term depends on ψ ( r )’s boundaryconditions. By choosing identical conditions and trivialvalues for the Green distribution function at the bound-aries ( G = 0 for Dirichlet or G ′ = 0 for Neumann) we arecapable of solving an infinite number of alike boundaryvalue problems.The discussion for the existence of the weight func-tion can be answered mathematically. Given the rela-tionship required, the weight is defined as w ( r ′ , r ) :=exp[ − γ ( r ′ , r )], yielding ~ ∇ { r ′ } γ ( r ′ , r ) = − ~f ( r ′ ). Consid-ering that the curl of the gradient of any scalar function istrivial then γ ( r ′ , r ) exists if and only if ~ ∇ × ~f = 0, whichmeans that ~f must be a conservative vector field! Withinthis view, γ ( r ′ , r ) represents the scalar potential associ-ated with a force. Anticipating this last restriction, thesolution for γ ( r ′ , r ) is independent of a path that simplyconnects r to r ′ yielding, w ( r ′ , r ) ≡ e − γ ( r ′ ) e − γ ( r ) = 1 w ( r , r ′ ) , (8)reflecting on the symmetry of the distribution as it willbe shown later. Finally, we are ready to define the inverseoperator of ˆ L asˆ L − { r } (cid:3) = Z r ′ G ( r ′ , r ) (cid:3) ( r ′ ) w ( r ′ , r )d r ′ , (9)conditioned by the boundary-value problem, which inturn defines G ( r ′ , r ) from eq. (4).There is one last piece of the puzzle to be resolved andit is related to the symmetry of the Green function distri-bution. Let us evaluate ˆ L { r } G ( r ′ , r )—i.e., the operatoracting on the second variable. Direct application of ˆ L { r } on eq. (7), using eq. (2), φ ( r ) = Z r ′ ˆ L { r } { G ( r ′ , r ) w ( r ′ , r ) } φ ( r ′ ) d r ′ + ˆ L { r } { b.c. } , (10)hints how this operator appears to work and leads us toanticipate the convolution of a Dirac distribution. Indeedthis is true. To clarify, here we exchanged integral andLiouville operators because they are acting on separatevariables, and the weight and Green functions (except at r ′ = r ) are differentiable.This conjecture can be proved from the following state-ment: two separate problems with different boundaryvalues and identical inhomogeneous differential equation—eq. (2)—share the same Green function distributionand satisfy eq. (10); therefore, by comparing equationsfor any two cases leads to ˆ L { r } { b.c. } = 0 because wecan always choose convenient trivial boundary values ( i.e. b.c. = 0) in one case. Consequently,ˆ L { r } { G ( r ′ , r ) w ( r ′ , r ) } ≡ δ ( r ′ − r ) , (11)an equality that bears meaning in the sense of the distri-butions. These final result unravels the symmetry of theGreen distribution function via the weight function, i.e. G ( r ′ , r ) = G ( r , r ′ ) w ( r , r ′ ) . (12)An interesting question now arises, and it is relatedto the possibility of using eq. (12) to drop the weightfunction out of the equation. This operation, with theaddition of the relation ~ ∇ { r ′ } w ( r , r ′ ) = − w ( r , r ′ ) ~f ( r ′ ),leads to, ψ ( r ) = Z r ′ G ( r , r ′ ) φ ( r ′ ) d r ′ + I ∂ r ′ h ψ ( r ′ ) ~ ∇ { r ′ } G ( r , r ′ ) − G ( r , r ′ ) (cid:2) ψ ( r ′ ) ~f ( r ′ ) + ~ ∇ { r ′ } ψ ( r ′ ) (cid:3)i · n d S ′ . (13)Notice that for Neumann boundary conditions (NBC),unlike Dirichlet (DBC), both ψ ( r ) and its derivative —atthe boundaries— are necessary. Ergo, eq. (13) is incon-venient for NBC unless either ψ vanishes or ~f ( r ) = 0. Insuch a case, it deems necessary to use the version thatincorporates weight function.Actually, the vector field ~f does not appear in someof the Liouville operators used in physics. For instance,the Green function associated with the electrostatic field satisfies the relation ~ ∇ G ( r , r ′ ) = − πδ ( r − r ′ )—the irrel-evant factor of − π appears by convenience. The staticregime of the Klein Gordon equation—which also leadsto the Yukawa potential—also follows a similar behavior,as its associated Green function in 2D is K ( µr ), satisfy-ing the DE ( ~ ∇ − µ ) G ( r ) = − πδ ( r ) [26]. Clearly, ~f isabsent in both systems.Yet the DEs describing the behavior of other physi-cal systems such as the driven damped harmonic oscilla-tor, the diffusion equation at thermal equilibrium withan anisotropic diffusion coefficient, and the electrostaticpotential in the presence of anisotropic media, includethe existence of a vector field ~f —see section IV for moredetails about the first system.Surprisingly, any dependence on the weight function ineq. (13) has vanished. As previously stated, the weightfunction can only be defined when ~f is a conservativefield. Then, an important question now arises: is eq. (13)still valid for non-conservative vector fields? This in afundamental question that can be addressed in a futurework. Since the main purpose is to present a compactand rigorous algorithm to solve the Green function in 2Dspace, we will restrict our analysis to only the supportedcases. B. Boundary conditions
As previously stated, the Green function convenientlyinherits identical types of conditions as the target func-tion ψ at the boundaries. These can be summarized as,DBC : → G ( r , r ′ ) | r at R ext / int = 0 , NBC : → ∂ r G ( r , r ′ ) | r at R ext / int = 0 . (14)However, there are two hidden additional conditions thatmust be satisfied enforced by the presence of Dirac’s dis-tribution. The rationale behind is that without them G = 0 will be a solution to the Green function for thesimple boundary value problem. While this is directlyvisible for Dirichlet, notice that it also applies for Neu-mann’s case. The added restrictions appear at the ar-tificial boundary r = r ′ implying continuity of G anddiscontinuity of the local derivative. Both are essentialto secure a non–zero solution. Continuity is often re-garded considering that the Green distribution is still afunction and its derivatives up to second order exist inthe classical sense of the DE everywhere except at r = r ′ .Though there is a stronger argument that stems from the The 1D driven damped harmonic oscillator is modeled by theDE m ¨ x + b ˙ x + kx = F ( t ). The damping constant b plays the roleof ~f in this one dimensional system. Although the time t is therelevant variable describing this system (instead of the position x ), the one dimensional formalism we describe is analog to thismodel. Figure 1: The contours used in path integration: to theleft S ǫ and S δ to the right. The radius and thickness arechosen purposely as δ < ǫ/ ǫ → r = r ′ , conditions are deriveddirectly from eq. (4) by integrating over r inside the vol-ume delimited by the surface S δ enclosing r ′ such thatit is contained inside a vecinity— V δ —of r ′ (see right offig. 1 for an artistic view). Using the divergence theorem,the condition simplifies to,lim δ → I S δ ~ ∇ { r } G ( r , r ′ ) · ˆ n d S = 1 , (15)where we have kept the leading contributing term whiletaking the limit. For the one dimensional case it yieldsthe relation G ′ ( r ′ > , r ′ ) − G ′ ( r ′ < , r ′ ) = 1. C. Connection with the Sturm–Liouville problem
The weight function, if existent, is able to transformthe Liouville operator into a self–adjoint differential op-erator. Notice that action of w ( r ′ , r ) on eq. (1), w ( r ′ , r ) ˆ L { r ′ } (cid:3) = ~ ∇ { r ′ } · (cid:2) w ( r ′ , r ) ~ ∇ { r ′ } (cid:3) (cid:3) + w ( r ′ , r ) g ( r ′ ) (cid:3) , (16)yields the otherwise known Sturm–Liouville form forPDE’s. Namely, the Sturm–Liouville differential oper-ator reads then as,ˆ L SL { r ′ } (cid:3) = w ( r ′ , r ) ˆ L { r ′ } (cid:3) . (17)Consequently, operating onto the Green distributionfunction gives equal results for both operators, i.e. ˆ L SL { r ′ } G ( r ′ , r ) = δ ( r − r ′ ).This last result connects the Sturm-Liouville problemwith null eigenvalues and the Green function distribu-tion problem where the former is the solution to the firststrictly when r = r ′ under either Dirithlet or Neumann boundary conditions. Resulting from this, G ( r ′ , r ) is con-tinuous everywhere and differentiable at r = r ′ ; the be-havior of its derivative at r = r ′ is dictaminated by theLiouville operator in the domain of the problem and spec-ified by the Dirac Delta distribution. D. Moving from 2D to 1D: the infinite coupling
As noted before, let us elaborate on the simplest sce-nario where a reduction of the dimension of the problemsignificantly improves our chances of procuring a gen-eral solution. Assume we would like to find the two-dimensional Green function in accordance with eq. (4).Taking advantage of the completeness of the Fourier in-finite expansion of any periodic function, we propose tosolve the two-dimensional DE in polar coordinates; thedimensional reduction occurs due to the periodicity in θ that does not take place in Cartesian coordinates.Although the Laplace operator in polar coordinates isknown to be separable in the variables r and θ , the in-troduction of the additional terms in eq. (1), as alreadymentioned, may lead to a DE that cannot be split con-veniently. Eq. (4) is then given by, h ∂ ∂r + 1 r ∂∂r + 1 r ∂ ∂θ + f r ( r, θ ) ∂∂r + f θ ( r, θ ) 1 r ∂∂θ + g ( r, θ ) i G ( r , r ′ ) = δ ( r − r ′ ) , (18)where convenient periodic conditions that must be sat-isfied suggests we should expand the Green functiondistribution in Fourier modes. Tentatively, we can re-sort to an expansion of the form P λ e iλ ( θ − θ ′ ) G λ tomatch the delta distribution expansion— i.e. δ ( r − r ′ ) = πr δ ( r − r ′ ) P λ e iλ ( θ − θ ′ ) —but since we cannot guaranteethat the G λ coefficients are θ ′ -independent (only under proper angular symmetry conditions) then we will as-sume G ( r , r ′ ) expands as, G ( r , r ′ ) = 12 π X λ ∈ Z e iλθ G λ ( r, r ′ , θ ′ ) . (19)With that in mind and multiplying eq. (18) by r to avoiddivergences at r = 0, X λ ∈ Z e iλθ h r d d r + r (1 + rf r ) dd r + ( − λ + iλr f θ + r g ) i G λ = X λ ∈ Z e i ( λ − λ ′ ) θ rδ ( r − r ′ ) . (20) The Fourier expansion, f ( r, θ ) = X µ ∈ Z f µ ( r ) e iµθ ; f µ ( r ) = 12 π Z π f ( r, θ ) e − iµθ dθ . Notice how we cannot obtain a solution because thereremains a residual dependence of θ in functions ~f and g . Despite this, a simplification can be manufacturedwhen they are replaced by their Fourier series form beforeintegrating on θ over a full period. This step yields ourmaster equation where we deduce that the G λ ( r, r ′ , θ ′ )modes satisfy the DE , r G ′′ λ + r G ′ λ − λ G λ + X µ ∈ Z r f r µ G ′ λ − µ + X µ ∈ Z (cid:2) ir ( λ − µ ) f θ µ + r g µ (cid:3) G λ − µ = rδ ( r − r ′ ) e − iλθ ′ . (21)This final result shows we have accomplished to reducethe rank of the effective Green function to solve at theexpense of requiring a countable large number of theseGreen modes. It is remarkable how the dependence on θ ′ is delegated to a quasi-negligible term at the right handside of the equation. We will develop this argument fur-ther in the following sections.This formulation represents an infinitely coupled sys-tem of linear differential equations that unsurprisinglycontains the solution to the Green function for the clas-sical source-free wave function; the structure of func-tions ~f and g defines the strength of the entanglement ofGreen’s free wave modes appearing in the rate at whichthe Fourier coefficients—functions—go to zero with in-creasing mode frequency. For simplicity, we opt to recallGreen’s Fourier modes as λ -modes, and ~f and g ’s modesas µ -modes suggested by the indexes employed in theequation above. E. More on boundary conditions of the λ -modes One last effort must be done to explain how boundaryconditions are inherited along the free-wave modes. Thekey to this understanding depends on the geometry ofthe problem and the originating expansion from eq. (19);we can identify two cases for disc–like systems: the an-nulus and the disc . Other geometries will be studied ina future work. For the annulus, either under Dirichlet orNeumann boundary conditions, the function or its deriva-tive must vanish at the boundaries. This can be met ifall modes preserve the vanishing values at both inner andouter boundaries —under uniform convergence. In doingso, we guarantee to meet all requisites for the Green func-tion and a solution is obtained. Conversely, preservingboundary conditions for the disc is not trivial becausewe do not have one but two boundaries (the second at r → + ). Due to the oscillating behavior of e iλθ with θ at r → + , all λ = 0 modes must vanish at the originto ensure continuity of the Green distribution function. We have dropped out the dependencies of all functions on r , θ , r ′ , and θ ′ facilitating a comprehensible reading. This can be enforced examining eq. (21) as r → + . Dis-continuity due to the source at r = r ′ may be neglectedfor now to realize that we can, while approaching the ori-gin, consider the behavior of each r , r , and r terms independently . We draw then conveniently, r (cid:2) − λ G λ (cid:3) = r (cid:1) ,r h G ′ λ + i X µ ∈ Z ( λ − µ ) f θµ G λ − µ i = r (cid:1) ,r h G ′′ λ + X µ ∈ Z f rµ G ′ λ − µ + X µ ∈ Z g µ G λ − µ i = r (cid:1) , , where we can choose, via r terms, that G λ = 0 ∀ λ = 0and G = 0. Plugging this sequentially into r and r terms hints G ′ λ = 0 and G ′′ λ = 0 ( ∀ λ = 0) assum-ing that lim r (cid:1) f θµ G λ − µ = 0 and lim r (cid:1) f rµ G ′ λ − µ =0 ∧ lim r (cid:1) g µ G λ − µ = 0 (with the exception of λ = 0where the µ = 0 term remains, thus we will choose G ′′ = − g G ) respectively.This is supported from continuity of g ( r ) everywhere inthe disc and from the definition of ~f ( r ), where limited bythe existence of w ( r ′ , r ), as the gradient of a scalar func-tion. If such function, γ ( r ), where to be free of patholo-gies and differentiable everywhere in the disc (includingthe origin) then lim r (cid:1) f rλ = 0 and lim r (cid:1) f θλ = 0 for λ = 0. It remains to say, that in order to fulfill all aboveconditions we will require that f θ and f r are finite as r → + . Looking under the hood of these assumptions,note that consequently the 0-mode has a logarithmic di-vergence when r ′ = 0, i.e. G ∝ r (cid:1) − log r .In summary, the conditions for the disc at the originare the following two only for r ′ > G λ (0 + , r ′ , θ ′ ) = 0 ∀ λ = 0 ,G ′ λ (0 + , r ′ , θ ′ ) = 0 ∀ λ . (22)Exceptions and particularities emerging from the specificform of functions ~f and g must be taken into accountwhen detailing the boundary conditions and may alterthe relationships obtained above.This relationship has to be completed with the result-ing relationship at the artificial boundary r = r ′ obtainedwhen using a complementary surface S ǫ correspondingto an open ring of ǫ > r ′ lim ǫ → Z π (cid:2) ∂ r G ( r , r ′ ) > − ∂ r G ( r , r ′ ) < (cid:3) d θ = 1 . (23)which entails the radial averaged contribution. We haveeliminated angular contributions by selecting the conve-nient contour S ǫ suggesting a pathway to extend it tothe λ -modes. Direct substitution of eq. (19) along witha convenient choice of unity—inspired by eq. (21)—givesus ultimately, r ′ (cid:2) G ′ λ ( r ′ > , r ′ , θ ′ ) − G ′ λ ( r ′ < , r ′ , θ ′ ) (cid:3) = e − iλθ ′ , (24)where primes denote partial derivatives with respect tothe first argument at both left ( < ) and right ( > ) handsides of r ′ . Substituting this relationship in the differ-ential equation, we obtain a similar relationship for thesecond derivatives, essential to the numerical method, asfollows,( r ′ ) (cid:2) G ′′ λ ( r ′ > , r ′ , θ ′ ) − G ′′ λ ( r ′ < , r ′ , θ ′ ) (cid:3) = − e − iλθ ′ [1 + r ′ f r ( r ′ , θ ′ )] . (25)For the disc, the r ′ = 0 case must be clarified. Inpolar coordinates, Dirac’s distribution is best describedas absent of angular dependence, which entails that forall λ -modes except λ = 0 it is exactly zero. Therefore,the boundary at the origin for each λ -mode is dictatedby symmetry except for λ = 0. This last, carries thelogarithmic divergence. This means that conditions fornon-zero modes are unchanged. For the zero mode anddue to symmetry G ′ | r → + = − G ′ | r → − and G ′′ | r → + = G ′′ | r → − ; however, due to the divergence a cutoff mustbe set in place. Such a choice of cutoff will be discussedlater. III. FINITE DIFFERENCES METHOD, FDMOR FEM ON A REGULAR GRID
The FDM, or uniform mesh FEM, has been used exten-sively in the literature to find approximate solutions formany physical systems and its stability makes it a suit-able candidate to obtain a numerical Green distributionfunction. Some examples include the one-dimensionalSchr¨odinger equation [28], the Poisson equation for Elec-trodynamics [29], the Euler equations of inviscid fluidflow [30], solutions to 1D and 2D Burgers’ equation [31]and the time-fractional diffusion equation [32]. From themathematical perspective, the same method has been im-plemented to solve elliptic, hyperbolic and parabolic par-tial DEs on irregular meshes [33], with interfaces [34], orin finding optimal algorithms on nontrivial meshes [35].Orchestrating an exact solution to eq. (21) is virtuallynot possible. There are four cases where an analytical ap-proach can be attempted: two cases where either ~f or g are zero, requiring to find a base of G λ ’s that can decou-ple the system—hence, a diagonalization—, the uniquecase where the same base applies to both coupling matri-ces accompanying G ′ λ and G λ , and the trivial free -wave( ~f and g zero). Excluding the latter, finding this diag-onalizing operator for the first three cases will be ad-dressed in a future study.Therefore, we will compute a numerical solution wherewe approximate the operator with finite differentiation(the finite difference method —FDM or FEM for a reg-ular grid) and bind expansions to include all relevantGreen and function modes up to a calculated cutoff; max-imum and minimum modes will be chosen respectively for λ - and µ -modes symmetrically as | λ | ≤ L and | µ | ≤ M Variable or function Equivalent array r and r ′ r j = r + h j with r = R int j ∈ { , , , . . . , N } G λ ( r, r ′ , θ ′ ) G j,kλ | θ ′ = G λ ( r j , r ′ k , θ ′ ) P ( r ) := r P j = P (cid:0) r j (cid:1) Q µ ( r ) := r f r µ ( r ) + rδ µ Q jµ = Q µ (cid:0) r j (cid:1) R λ,µ ( r ) := r g µ ( r ) − λ δ µ R jλ,µ = R λ,µ (cid:0) r j (cid:1) + ir ( λ − µ ) f θ µ ( r ) Table I: A summary on the change of notation fromcontinuous to discrete formconsidering that M ≤ L for reasons that will be clarifiedafterwards.Since the Green function is twice differentiable, when r = r ′ , its Fourier series converges uniformly and its co-efficients decay at least as λ − , conditioned by equallywell behaved functions ~f and g . Then, a possible edu-cated choice of L is the minimum integer such that thesum of 1 /k up to L exceeds π p , with p a percentage ofaccuracy; for example, to achieve at most 1% of estima-tion error we require L >
A. A large matrix equation
The Finite Differences Method (FDM or FEM—finiteelements method—with uniform grid) is a simple ap-proach to computing derivatives of functions at a pointby using Taylor expansions on a discretized mesh. Indoing so, a derivative will rely on knowledge of the val-ues of the function in neighboring sites. Such is the art ofcomputing derivatives. The number of neighboring sitesto be taken into consideration determines the degree ofwhich the function approaches to the point value. Forinstance, in the so called three-point stencil (the site inquestion and its two adjacent neighbors), the first andsecond derivatives are accurate up to order square of themesh size. The choice of whether dissecting uniformly or non-uniformly ishighly dependable on the problem. For example, if we were in-terested in fracture dynamics we would prefer a non-uniform gridto model complex material topologies.
Going back to our problem in eq. (21), we turn toa simply redefined one dimensional DE for a sketch ofthe forthcoming operations. The left hand side readsrewritten as, P ( r ) d G λ ( r, r ′ , θ ′ )d r + X µ ∈ Z Q µ ( r ) d G λ − µ ( r, r ′ , θ ′ )d r + X µ ∈ Z R λ,µ ( r ) G λ − µ ( r, r ′ , θ ′ ) . In transforming the continuous variables r and r ′ intoa discrete equally–spaced mesh of size h , we will adoptmatrix notation for variables and functions; ergo, for N partitions defining N + 1 points h = ( R ext − R int ) /N ina disc–like geometry. For clarity, we summarize nota-tion changes in table I. This procedure applied over the aforementioned equation gives for r = r ′ , P j (cid:20) h X η ∈A | j a (2) η | j G j + η,kλ | θ ′ (cid:21) + X µ ∈ Z Q jµ (cid:20) h X η ∈A | j a (1) η | j G j + η,kλ − µ | θ ′ (cid:21) + X µ ∈ Z R jλ,µ G j,kλ − µ | θ ′ + O ( h ξ ) , with ξ the order of approximation, A | j the set of neighborsite indices, and a ( n ) η | j the respective coefficient (namelythe finite difference coefficient included into a matrix rep-resentation A ( n ) —see section B) of the η -th neighbor re-quired to compute the n -th derivative up to a predeter-mined order of accuracy [36]; in the three-point stencilcase, ξ = 2. Both sets of neighbor indices and coefficientsdepend on the information of the site j under inspection;if for example we are at or near an interface, boundaryor discontinuity then the strategy for choosing neighborsmay differ; we might be interested in computing deriva-tives using only points in regions where it makes sense.With some reorganization, the generated discrete DEcan be regarded as a matrix multiplication. To seethis, first we realize that by understanding G j,kλ | θ ′ as the( j, k )-th element of a constructed matrix G λ | θ ′ —of size N + 1 × N + 1—we can envision a column matrix vector G θ ′ that contains all λ -modes, or all of { G λ | θ ′ ∀ λ ∈ Z } ,where all operations from the previous complex arrayequation are condensed into an equally conceived matrix U multiplying G θ ′ . The following is a view of G θ ′ , G θ ′ = ... G − L | θ ′ ... G − | θ ′ G | θ ′ G | θ ′ ... G L | θ ′ ... , with G λ | θ ′ = G , λ | θ ′ G , λ | θ ′ · · · G ,N − λ | θ ′ G ,Nλ | θ ′ G , λ | θ ′ G , λ | θ ′ · · · G ,N − λ | θ ′ G ,Nλ | θ ′ ... ... . . . ... ...G N − , λ | θ ′ G N − , λ | θ ′ · · · G N − ,N − λ | θ ′ G N − ,Nλ | θ ′ G N, λ | θ ′ G N, λ | θ ′ · · · G N,N − λ | θ ′ G N,Nλ | θ ′ . (26)In principle, both matrices are infinitely large but forpractical terms they will be truncated on both λ - and µ -modes as mentioned in the previous section. Despitethis numerical simplification that will be carried out inthe numerical analysis, the infinite matrix U has a welldefined structure as will be detailed in section III B.Finally, the terms to the right of eq. (21) vanish for all r = r ′ leading us to believe that if U is invertible then thesolution to the discrete Green function G θ ′ is identicallyzero. However, attention should be paid at r = r ′ for its effect discards the trivial solution. Along with the othergeometrical boundary conditions the problem will nowhave a unique solution. These boundary conditions willbe addressed in section III C. B. Infinite matrix U To understand the structure of U we turn to the setof operations for a particular λ -mode. Seeing as U isinfinite we may encode rows by the integer value of the mode being solved and columns by the value of the modebeing correlated. Thus taking row λ from U , U λ G θ ′ = · · · , column λ − µ z }| { h Q µ + R λ,µ , · · · , ··· , M , ··· , , , ←− µ | column λ z }| { h P + 1 h Q + R λ, µ =0 , · · · , column λ − µ z }| { h Q µ + R λ,µ , · · · | µ −→ − , − , ··· , − M , ··· ... G λ − µ | θ ′ o row λ − µ ... G λ | θ ′ o row λ ... G λ − µ | θ ′ o row λ − µ ... , (27)with the following definitions for matrices P , Q µ , R λ,µ , P = P · · · P · · · ... ... . . . ... · · · P N × A (2) , (28) Q µ = Q µ · · · Q µ · · · ... ... . . . ... · · · Q Nµ × A (1) , (29) R λ,µ = R λ,µ · · · R λ,µ · · · ... ... . . . ... · · · R Nλ,µ . (30)One last remark on matrix U is that the density of non–zero entries is at most 3 /N for the three-point stencil.For N sufficiently large, it will become essential to finda way to manage such sparsity for all speedups, data-compression and efficiency in memory footprint. C. Discrete Boundary Conditions
Retaking conditions detailed thoroughly in section II Band at end of section II D we are now in capacity of pa-rameterizing the values of G j,kλ | θ ′ . This parametrizationshould further reflect the behavior of the δ –function. Thefollowing are the conditions for the 3–point–stencil: ( i. )at r = R ext , G N,kλ | θ ′ = 0 for DBC ; G N +1 ,kλ | θ ′ − G N − ,kλ | θ ′ = 0 for NBC , (31) ( ii. ) at R int for the annulus, G ,kλ | θ ′ = 0 for DBC ; G ,kλ | θ ′ − G − ,kλ | θ ′ = 0 for NBC , (32)( iii. ) for the disc at R int (disregarding G ′ λ = 0 for now), G ,k | θ ′ − G − ,k | θ ′ = 0 λ = 0 ,G ,kλ | θ ′ = 0 λ = 0 , (33)and, finally, ( iv. ) at the interface r = r ′ the conditionreads, (cid:0) G k > +1 ,kλ | θ ′ − G k > − ,kλ | θ ′ (cid:1) − (cid:0) G k < +1 ,kλ | θ ′ − G k < − ,kλ | θ ′ (cid:1) = 2 hr k e − iλθ ′ . (34) Here we have adopted the subscript convention of <, > torefer to points to the left and right of the site of derivativeevaluation. Note how all equations above reference andhighlight a few fictitious points. The mesh points thatlay outside or beyond the valid grid are G N +1 ,kλ | θ ′ , G − ,kλ | θ ′ , G k < +1 ,kλ | θ ′ , and G k > − ,kλ | θ ′ . These spurious terms must bedealt with and simplified in order to be able to incorpo-rate readily all conditions. D. A Non–Trivial Matrix Equation and a Solution
We will now show the explicit matrix equation associ-ated with the conditions described above. As mentioned,they depend on the degree of accuracy that we choose,or equivalently, the stencil. We will describe the proce-dure for the three-point stencil and further discuss howto generalize for higher orders of approximation.With the boundary relationships in mind, here ineqs. (31) to (34), eq. (21) (multiplied by − h ) equatespartially to zero (when r = r ′ ) as,2 P j G j,kλ | θ ′ − P j ( G j +1 ,kλ | θ ′ + G j − ,kλ | θ ′ ) − h X µ R jλ,µ G j,kλ − µ | θ ′ − h X µ Q jµ ( G j +1 ,kλ − µ | θ ′ − G j − ,kλ − µ | θ ′ ) = 0 , where via eq. (34) the latter can be used to simplifyboth spurious terms (appearing at r = r ′ ) G k < +1 ,kλ | θ ′ and G k > − ,kλ | θ ′ . After crossing out these terms by iterative sub-stitution we obtain a generalized expression for the abovevalid for almost every point in the grid. The general dis-crete equation yields for r ′ > P j G j,kλ | θ ′ − P j ( G j +1 ,kλ | θ ′ + G j − ,kλ | θ ′ ) − h X µ R jλ,µ G j,kλ − µ | θ ′ − h X µ Q jµ ( G j +1 ,kλ − µ | θ ′ − G j − ,kλ − µ | θ ′ ) = − hr j δ j,k e − iλθ ′ × (cid:26) − h r k ) [1 + r k f r ( r k , θ ′ )] (cid:2) − δ λ, (cid:3)(cid:27) , (35)where the new term that accounts for the boundary con-dition at r = r ′ has appeared. Due to the absence ofa left-hand limit as r ′ = 0, according to eq. (34), thisterm is exactly − hr j δ j,k e − iλθ ′ at the origin. Actually,this condition holds for the mode λ = 0 in general dueto translational invariance—this invariance is clearly ab-sent for the other modes. The terms composing the righthand side of last equation can be viewed as of order ofmesh–size or order of radial distance from the origin asfollows,1. − h f r ( r k , θ ′ ), a surprising third order correctiondue to the vector field appearing after substitutingthe interface difference in derivatives.2. h r k = h R int + h k , the leading order that sub-stituted yields a first order constant term and asecond order increasing term.3. − h r k = − h R int + h k , a negative term significantcloser to the origin. As expected, the behavior ofthe discrete version near zero validates our previouschoice of boundary condition for the disc.4. For r ′ = 0, we must implement a cutoff such that r = ǫ > ǫ will be discussed below.Because the error in the differential equation is of O ( h ),we should incorporate all terms to the calculation. How-ever, we will neglect the higher order term—first term—since this will simplify our calculations of ψ ( r ).This final expression is valid everywhere including thecontroversial j = 0 , N points, where either Dirichlet orNeumann conditions complete eq. (35) at the borders. Inthose two cases, substitutions must take place following Mode U j,jλ,λ U j,j ± λ,λ U j,jλ,λ − µ U j,j ± λ,λ − µ ∀ λ P j − h R jλ, − P j ∓ h Q j − h R jλ,µ ∓ h Q jµ Table II: Nonvanishing matrix elements of U for1 ≤ j ≤ N − U corresponding to exterior and interior bor-ders are modified. See the substitution rules in tables IIand III.We now define our complete matrix system as U · G θ ′ = − h V · E θ ′ . Among other things, the right hand side ac-counts for the contribution of Dirac’s distribution. Thetwo additional definitions appearing correspond to firsta distance parameter generalized into α jλ , a new objectthat incorporates the boundary conditions at both j = 0and j = N . Notice, for instance, that keeping the term r j at every point does not explain the vanishing of theGreen function at the boundaries when DBC are con-sidered, neither does it describe the correct behavior at r = 0 for a disk. Actually, when the last condition is con-sidered, an ultraviolet cutoff ǫ —such that ǫ → r = r ′ and a logarithmic behavior near the origin when r ′ = 0.Although the appearance of this divergence can easilybe visualized after studying the behaviour of eq. (35) at j = 0 for the mode λ = 0 in a disk, its existence at anypoint—also for an annulus—is guaranteed by the infinitenumber of λ -modes that must be summed up to obtainan exact solution. Therefore, it is not surprising that ǫ and h are related—see section IV for more details.Matrix terms are written as, V j,kλ,µ = α jλ δ λ,µ δ j,k , (36a) E j,kλ | θ ′ = e − iλθ ′ δ j,k , (36b)and the solution to the λ –modes matrix G θ ′ is, G θ ′ = − h A · E θ ′ , (37)where we have defined A = U − · V assuming that U isinvertible.Matrix A was declared because it has interesting sym-metry properties that will be discussed in the next sec-tion. Tables II to V outline how to fill the matrix ele-ments of the objects we have described.0 Mode j U j,jλ,λ U j,j +1 λ,λ U j,j − λ,λ U j,jλ,λ − µ (D) ∀ λ , N ∀ λ P − h R λ, − P N/A − h R λ,µ (N) ∀ λ N P N − h R Nλ, N/A − P N − − h R Nλ,µ
Table III: Nonvanishing matrix elements of matrix U that define DBC (D) and NBC (N) for an annulus.N/A specifies those elements that lay outside U . Mode j U j,jλ,λ U j,j +1 λ,λ U j,j − λ,λ U j,jλ,λ − µ λ = 0 0 2 − h g − λ = 0 N λ = 0 0 , N λ = 0 0 2 − h g − ∀ λ N P N − h R Nλ, N/A − P N − − h R Nλ,µ λ = 0 0 1 0 N/A 0 Table IV: Nonvanishing matrix elements of matrix U that define DBC (shown above) and NBC (shownbelow) for a disc. N/A specifies those elements that layoutside U . E. The parameter α and the symmetry of A A closed relation can be found for the matrix describ-ing the entire Green function. Using the results fromprevious section it is G j,kθ,θ ′ = − h π X λ, µ e iλθ e − iµθ ′ A j,kλ,µ , (38)where A j,kλ,µ = [ U − ] j,kλ,µ α kµ . As previously mentioned, theparameter α kλ generalizes the radial parameter r k , includ-ing the boundary conditions. On the other hand, it isworthwhile to state the symmetry conditions that A sat- Mode j α jλ (DBC) α jλ (NBC) Geom. λ = 0 0 ǫ − ǫ − (D) λ = 0 0 0 r (A) λ = 0 1 ≤ j ≤ N − r j r j (A, D) λ = 0 N r N (A, D) λ = 0 0 0 0 (D) λ = 0 0 0 r − h r (A) λ = 0 1 ≤ j ≤ N − r j − h r j r j − h r j (A, D) λ = 0 N r N − h r N (A, D) Table V: Elements α jλ that define DBC and NBC for anannulus (A) and a disc (D). isfies Im (cid:16) A j,k , (cid:17) = 0 , A j,k − λ, = A j,kλ, , A j,k , − µ = A j,k ,µ , A j,k − λ,µ = A j,kλ, − µ , A j,k − λ, − µ = A j,kλ,µ . (39)Using table V, it is easy to see that the matrix elements[ U − ] j,kλ,µ satisfy the same symmetry properties. F. The algorithm
The algorithm for a numerical solution can be summa-rized as follows:1. The values of L and M are determined accordingto the required level of approximation.2. We fill all elements described in table I; the matrixelements U j,kλ,µ are filled by blocks using the rulesshown in tables II, III, IV. Matrix elements α jλ arealso filled according to table V.3. Matrix U is inverted and so matrix A is computed.4. The Green function is computed according toeq. (38).Using previous results, we can deduce a closed form for ψ ( r ) for both DBC and NBC using the conventions statedin eq. (7) and eq. (13). Although there are many ways toperform the integrals stated in previous equations, andthe reader can choose the method that he or she prefers,a sketch of these solutions, using the trapezoid rule, isshown in appendices F and G. Particular cases and properties
We will analyze some particular cases that can be de-duced from the procedure explained above. We will startfocusing on the one-dimensional case.
One dimensional case
The analysis of a Green function in one dimension re-quires an appropriate definition of a general DE obeyedby the Green function G ( x, x ′ ). Unfortunately, a directanalysis of the results by studying eq. (21) is not straight-forward due to the clear differences between the Lapla-cians in cartesian and polar coordinates. Let us imaginea general second order DE of the form L x ψ ( x ) = φ ( x ),where the Green function satisfies the relation L x G ( x, x ′ ) = P ( x ) d Gdx + Q ( x ) dGdx + R ( x ) G = βδ ( x − x ′ ) . (40) A more detailed derivation can be found in appendix D; z denotesthe complex conjugate of z . P ( x ) might not be necessary, as it can beeliminated by division, but its inclusion allows us to havea more general analysis. The constant therm β seemsclumsily placed, as its value is usually 1. Nonetheless,some formalisms define the Green function by means ofthe operator L x G ( x, x ′ ) = − δ ( x − x ′ ), thus introducinga change of sign that can be contemplated in our study.By following a similar analysis as that shown above,we can deduce an appropriate recurrence relation foreq. (40), which is(2 P j − h R j ) G j,k − ( P j + h Q j ) G j +1 ,k − ( P j − h Q j ) G j − ,k = − hβ δ j,k . (41)Having confined the system within the domain x ∈ [ x , x N ], our step size is now h = ( x N − x ) /N .From this point on, we can apply the results obtainedfor the two dimensional problem in this study. Noticethat, in the absence of modes that account for the angulardependence, we can always say that A j,kλ,µ = A j,kλ,µ δ λ δ µ .Therefore, eq. (38) becomes G jk = − h A j,k , where A j,k = [ U − ] j,k α k (42)and α k accounts for the boundary conditions. By makingthe association x j = x + hj , the elements described ineq. (42), for 1 ≤ j ≤ N −
1, are now filled using thefollowing rules:1. U j,j = 2 P j − h R j .2. U j,j ± = − P j ∓ h Q j .3. U , = U N,N = 1 for DBC. For NBC: U , =2 q − h b , U N,N = 2 P N − h R N , U , = − P and U N,N − = − P N − .4. α j = β .5. α = α N = 0 for DBC. For NBC: α = β and α N = β .The weight function, which now guarantees the sym-metry condition G kj = w jk G jk takes the form w ( x ′ , x ) = e U ( x ′ ) e U ( x ) , U ( z ) = 1 P ( z ) exp (cid:20) Z zz Q ( y ) P ( y ) dy (cid:21) , (43)where z is an irrelevant constant. Solutions for ψ ( x )with both DBC and NBC using the trapezoid rule asmethod of integration are shown in appendix G. Monopole–like case
This takes place when both ~f ( r ) and g ( r ) have no sig-nificant angular dependence, so the mode µ = 0 is theironly relevant contribution; this implies that Q jµ = R jλ,µ =0 for µ = 0. The off-diagonal matrices U j,kλ,λ − µ thusvanish—this leads to a block diagonal U matrix—and so the system becomes separable in the radial and angu-lar variables. Having now the relation A j,kλ,µ = A j,kλ,µ δ λ,µ ,eq. (38) reduces to G jkθθ ′ = − h π X λ e iλ ( θ − θ ′ ) A j,kλ,λ . (44)Notice that each mode can now be solved independently. G. Beyond the three-point stencil
As mentioned, we only showed an explicit analysis fora three-point stencil approximation. This method canbe generalized to include the contribution of more neigh-bors in the derivative terms, i.e., higher order stencilsthat provide more accurate degrees of approximation in h . In spite of its simplicity, the three-point stencil has thegreat advantage that the spurious terms that arise fromthe boundary conditions can be eliminated in a simplefashion.The description of the system with a five-point stencil,for instance, will increase the amount of terms differentfrom zero in U —for example, terms of form U j,j ± λ,λ − µ willprovide non-zero contributions. Having a higher degreeof approximation, that demands the inclusion of morenon-trivial matrix terms, the grid size can be reduced.Although there is no guarantee that the inversion pro-cess is optimized in time when the contributions of moreneighbors are included, as the matrices are highly sparse,there is a clear optimization of memory storage.Yet a great disadvantage that higher stencils inheritis the elimination of the spurious terms that come fromthe boundary conditions. For instance, when we dealwith the condition at j = k ( r = r ′ ), eq. (34) will in-clude more coefficients outside the grid, so the recur-rence relation that is obtained will not be able to elimi-nate all of them—at least, using the same procedure weimplemented. Therefore, a different approach must beperformed. A possible solution could be expanding thederivatives around a point different from the center, soavoiding the spurious terms; this process is studied indetail in [36]. Nonetheless, this could be discussed in afuture study. IV. NUMERICAL RESULTS
We will use the formalism described above to solvesome particular examples.
Example 1: A one dimensional case
As a first example, let us study a one dimensional sys-tem with a known analytical solution, useful to test theformalism we have described. Let us suppose we want tosolve the DE in the domain [0 , x ψ ′′ ( x ) + xψ ′ ( x ) + ( x − ψ ( x ) = J ( x ) , (45)2 Func. ψ N =32 ψ N =128 ψ N =512 ψ DBC max . . . . × − . × − . × − MSE 3 . . × − . × − ψ DBC max . . . . . × − . × − MSE 3 . . × − . × − f NBC min − . − . − . . × − . × − . × − MSE 1 . × − . × − . × − Table VI: Analysis of the accuracy of the two numericalsolutions to eq. (45). The first two blocks use theconditions ψ (0) = 0 and ψ (4) = 2 and show the value ofthe maximum—exact value ψ DBC max = 2 . . . . , itspercentage error (PE), and mean square error (MSE) ofthe function along the domain [0 , ψ ′ (0) = 0 and ψ ′ (4) = 4 and focus on the globalminimum —exact value ψ NBC min = − . . . . usingthe weight function—remember that for NBC theformalism with the weight function is required.where J n ( x ) and Y n ( x ) are the Bessel functions of firstand second kind of order n . Using eq. (43), we can easilydeduce that w ( x ′ , x ) = x/x ′ .The conditions ψ = 0 and ψ N = 2 (Dirichlet), lead tothe analytical solution, ψ ( x ) DBC = 124 J (4) [ J ( x )(48 − J (4)++ πxJ (4) J ( x ) Y ( x )) − πxJ (4) J ( x ) J ( x ) Y ( x )] . (46)Conversely, with conditions f ′ = 0 and f ′ N = 2 (Neu-mann), the analytical solution yields, ψ ( x ) NBC = 124 x ( J (4) − J (4)) [(2 J (4) + 3 J (4)) × ( x ( x − J ( x ) − x − J ( x ))+ x J ( x )(96 + 2 J (4) − J (4))] . (47)Analytic and numerical results are compared for bothcases in fig. 2 and table VI—see eq. (G4a) and eq. (G4b)for explicit expressions using a numerical approach.It is interesting to contrast the numerical solutionsshown in fig. 2 using the weight function and the onethat arises without the weight function formalism—performing the replacement w jk G jk → G kj . Interest-ingly, from table VI we conclude that the introduction ofthe weight function leads to a more accurate result. − − − − − − − − Figure 2: Above: solution to the DE given by eq. (45)with the initial conditions ψ (0) = 0 and ψ (4) = 2.Below: solution to the DE given by eq. (45) with theinitial conditions ψ ′ (0) = 0 and ψ ′ (4) = 2. Table VIanalyzes the accuracy of the numerical solutions. Note:we made x = 10 − to avoid numerical divergences at x = 0. Example 2: The two-dimensional Helmholtzequation with imaginary wave number
We now shift our attention to solve a two dimensionalsystem. Let us consider the DE( ~ ∇ − m ) G ( r , r ′ ) = δ ( r − r ′ ) . (48)In the absence of the vector field ~f ( r ), we conclude that w ( r , r ′ ) = 1; besides, the system is separable in the radialand angular coordinates. The solution to last equationconfined in a large disc of radius r N = R ext with Dirichletboundary conditions can be found analytically. Adapt-ing the result found in [13], we deduce that the Greenfunction associated with eq. (48) is G ( r , r ′ ) = − m π X λ e iλ ( θ − θ ′ ) h I λ ( mr < ) K λ ( mr > ) − t λ ( R ) I λ ( mr ) I λ ( mr ′ ) i , (49)where r > and r < are the maximum and minimum be-tween r and r ′ , I λ ( x ) and K λ ( x ) are the well-known3 (cid:1) (cid:0) (cid:2) (cid:3) (cid:4) (cid:5) Figure 3: Solution to eq. (48) using the numericalsolution explained in section III (continuous gray lines)and using eq. (49) (black dashed lines) for θ = θ ′ anddifferent values of r ′ . We have chosen in all cases unitssuch that m = 1. We also used L = 80 for both cases.In the numerical solutions N = 4096 and ǫ = 0 . h , inthe analytical solution s = 0 . h . The values of theminima and the Mean square error (MSE) are shown intable VII.modified Bessel function of the second kind and t λ ( R ) = K λ ( mR ) I λ ( mR ) . Taking a look to eq. (49) we deduce thatthe Green function diverges—many distributions are for-mally infinite. Actually, the first term of previous sumcan be reduced to − m π K ( m | r − r ′ | ).Notice that the Green function diverges logarithmi-cally as as r = r ′ , as K ( x ) x → ∼ ln(2 /x ) − γ , with γ the Euler Mascheroni constant. A cutoff s , which rep-resents a minimum separation distance between r and r ′ [11], is usually introduced to address this divergence. Inturn, s might be related to L and ǫ .We now implement the numerical analysis to verify lastsolution noticing that P j = ( r j ) , Q jµ ( r ) = α j = r j δ µ, ,and R jλ,µ = − ( λ + m P j ) δ µ, .The following step is determining an appropriate valuefor ǫ . From fig. 3 we see that the solution for a fixedvalue of r ′ shows a peak at r = r ′ . As we sum up allmodes the magnitude of the height of the peaks must beinfinite. However, the introduction of L guarantees thepeaks to be finite. The height of the peak at r = r ′ = 0is associated with ǫ , as the functions K λ ( r ) diverge at r = 0. The cutoff ǫ is chosen in such way that the heightof the peaks in the neighborhood of r = r ′ = 0—i.e., | r − r ′ | = O ( h )—are close enough. We found that for N ≥ , ǫ ≃ . h . Surprisingly, we found that ǫ doesnot depend on L for large enough L .We recall that this is an approximation; the exact so-lution for the distribution is found in the limits ǫ → L → ∞ . Comparisons between the numerical andanalytical results are shown in fig. 3 and table VII. Values of the minima G ( r, r ′ ) r ′ = 0 r ′ = 1 . r ′ = 3 . r ′ = 7 . − . − . − . − . − . − . − . − . . . ( − . ( − . ( − MSE 3 . ( − . ( − . ( − . ( − Table VII: Values of the minima shown in fig. 3 for theNumerical solution (NS) and analytical solution (AS).PE means percentage error (from 0 to 100%) and MSEis the mean square error. x ( y ) stands for x × y . Figure 4: Solution to ψ ( a ) ( r, θ ) for different values of θ in units in which m = 1. We used N = 256.Last result will now be used to solve a inhomogeneousequation of the form ( ~ ∇ − m ) ψ ( r ) = φ ( r ), whose gen-eral solution is provided in appendix G with r = R int =0. Now, let us consider φ ( r ) to be a function defined overa disc of radius R = 10 and study the two following cases:(a) ψ ( a ) ( r, θ ), with φ ( r ) = − r sin θ and ψ ( R, θ ) = 2.(b) ψ ( b ) ( r, θ ), with φ ( r ) = ( r − / , ≤ θ < π − r − / , π ≤ θ < π and ψ ( R, θ ) = ( , ≤ θ < π − , π ≤ θ < π .Solutions to ψ ( a ) ( r, θ ) and ψ ( b ) ( r, θ ) for some angles areshown in fig. 4 and fig. 5, respectively. Example 3: A pedagogical example
Now let us apply the same formalism to solve an-other two-dimensional problem. Let us suppose that4
Figure 5: Solution to ψ ( b ) ( r, θ ) for different values of θ in units in which m = 1. The continuous line shows thesolution for θ = { π , π , π , π } , the dotted line for θ = {− π , − π , − π , − π } . We used N = 256 and L = 80. Figure 6: Solution to G ( r , r ′ ) with DBC, as given inexample 3 for different values. We set h = andmade L = 15.we want to find the Green function associated with thetwo-dimensional DE [ ~ ∇ + ~ ∇ Z ( r )] · ~ ∇ ψ ( r ) = φ ( r ), where Z ( x, y ) = 2 x y . After transforming the system to po-lar coordinates, we can see that f r = r (1 − cos 4 θ ) and f θ = r sin 4 θ . The elements defined in table I now be-come Q jµ = r j δ ,µ + ( r j ) (cid:2) δ ,µ −
12 ( δ ,µ + δ − ,µ ) (cid:3) , (50a) R jλ,µ = − λ δ ,µ + 12 ( r j ) ( λ − µ )( δ ,µ − δ − ,µ ) . (50b)The system will be confined in an annulus or internalradius R int = 1 and external radius R ext = 2. Fig. 6and fig. 7 show the Green function for DBC and someparticular parameters but different values of L . Fig. 8shows the results for different parameters under NBC. (cid:6) (cid:7) (cid:8) Figure 7: Solution to G ( r , r ′ ) with DBC, as given inexample 3 for different values. We set h = andmade L = 30. (cid:9) (cid:10) (cid:11) (cid:12) (cid:13) (cid:14)(cid:15) (cid:16) Figure 8: Solution to G ( r , r ′ ) with NBC, as given inexample 3 for different values. We set h = andmade L = 40.Notice from fig. 6 and fig. 7 how the Green functions be-come zero at the borders and present a discontinuity at r = r ′ . As expected, the larger L , the larger the mag-nitude of the value at that point; however, this valuedecreases as | θ − θ ′ | increases. For NBC, as illustrated infig. 8, something similar happens. However, the deriva-tives at the borders are now the ones which tend no bezero. While this is clear as r →
2, the asymptotic behav-ior toward zero close to the inner border can be appreci-ated.When using the Green function to solve a particularinhomogeneous equation, it is clear that ~ ∇ × f = 0, sothe weight function exists. It is w ( r ′ , r ) = e r ′ (1 − cos 4 θ ′ ) e r (1 − cos 4 θ ) = P λ w λ ( r ′ ) e iλθ ′ w ( r ) , (51)whose only nonvanishing modes in discrete coordinates5are given by w j λ = ( − λ e ( r j ) I λ ( ( r j ) ). Although weare not interested in finding ψ ( r ) for a particular bound-ary problem, all the steps are carried out to accomplishthis goal. Example 4: The Stationary Diffusion Equation
It is worthwhile to describe how our formalism canbe adapted to solve the diffusion equation at “thermal”equilibrium. Let ψ ( r ) and D ( r ) represent the density ofthe diffusion material and the anisotropic diffusion coef-ficient, respectively. In the stationary regime, ψ ( r ) satis-fies the DE D ( r ) ~ ∇ ψ ( r ) + ~ ∇ D ( r ) · ~ ∇ ψ ( r ) = 0 . (52)Although eq. (52) does not have the standard form shownin eq. (2), after dividing eq. (52) by D ( r ) and defining ~f ( r ) as ~f ( r ) = D ( r ) ~ ∇ D ( r ), the standard form can beobtained. The weight function is now guaranteed toexist, as ~ ∇ × ~f = − D ~ ∇ D × ~ ∇ D + D ~ ∇ × ~ ∇ D = 0.Actually, w ( r , r ′ ) = D ( r ) D ( r ′ ) .The viability of our method to solve eq. (52) dependson the particular form of the diffusion coefficient. Thefollowing possibilities may arise: (a) ~f has no poles in thetwo dimensional domain; (b) ~f has a divergence in r = 0that can be eliminated once ~f is multiplied by r ; (c) thedivergence at r = 0—or any other radial divergence—previously discussed still persists after multiplication by r ; and (d) D − has poles for some θ ∈ [0 , π ).The cases (a) and (b) can be solved with the regularprocedure we have described; the modes f r µ and f θ µ arewell–behaved and so the elements Q jµ and R jλ,µ definedin table I exist. The possibility stated in (c) demands aredefinition of ~f to eliminate any possible radial diver-gence; however, this redefinition does not guarantee theexistence of the weight function. The situation describedin (d) is problematic, as some of the modes f r µ and f θ µ are divergent. Last situation is alleviated by workingwith the original DE, eq. (52); nonetheless, the processthat we must follow to solve a system whose mathemat-ical form differs from eq. (2) has not been described inthis work.Similar analysis can be performed as we deal with thePoisson’s equation associated with electrostatic potentialin an anisotropic media, among others. V. CONCLUSIONS
In this paper we have analyzed the Green function for-malism and studied under which conditions such mecha- nism can be used to obtain the solution of an inhomoge-neous DE. Particularly, we found that there exists a func-tion, which we called the weight function, that makes theLiouville operator self-adjoint. This function also definesthe symmetry properties of the Green function (how it istransformed under the exchange of r and r ′ ).After decomposing the Green function as a sum ofFourier modes, an infinite set of coupled second order dif-ferential for the radial variable is found. While such setdecouples when the initial DE is separable in the radialand angular variables, the coupling in the modes arisesas the vector field f ( r ) and the scalar function g ( r ) areexpressed as a sum of Fourier modes.An algorithm to solve the Green function associatedwith a general class of Liouville operator was solved usinga FEM. We used a simple three-point stencil approach toapproximate the solution and focused on both Dirichletand Neumann boundary conditions. A set of approxi-mations was made, which included a truncation of theinfinite number of modes, a minimum distance when thesystem is confined in a disc, and the discard of the term − h f r ( r k , θ ′ ). While the first two approximations arewell–justified because the Green function has a naturaldivergence, the last one was performed by convenience(anyway, it provides a very small contribution).The algorithm was verified by comparing with knownresults and obtaining very small percentage errors. Anadditional example whose solution cannot be found bymeans of the regular algorithms was shown.We consider that the presented method is a useful at-tempt to solve Green functions of operators whose ra-dial and angular variables cannot be separated. How-ever, we expect this algorithm can be improved by otherauthors in the future to obtain more accuracy withoutthe need of creating huge matrix systems, which demandlarge storage memory and computational time. Some ofthe improvements may include the implementation of themethod for higher order stencils, an optimized calcula-tion either mathematically or numerically of ǫ , simplifiedformulas for the calculation of ψ ( r ) or “on the go” algo-rithms that do not require the inversion of the matrix orthe storage of temporal information. ACKNOWLEDGMENTS
This work was partially funded by UniversidadCat´olica de Colombia. The diffusion coefficient is assumed to be well–behaved withinthe annulus or disc. However, ~f can have poles within the same domain. [1] J. Schwinger, Journal of Mathematical Physics , 407 (1961),https://doi.org/10.1063/1.1703727.[2] J.-S. Wang, B. K. Agarwalla, H. Li, and J. Thingna, Frontiers ofPhysics , 673 (2014).[3] S. Foster and N. Neophytou, Computational Materials Science , 91 (2019).[4] L. P. Kadanoff and G. Baym, “Quantum statistical mechanics.Green’s function methods in equilibrium and nonequilibrium prob-lems.” New York: W. A. Benhamin, Inc. XI, 203 p. (1962). (1962).[5] R. Alkofer and L. von Smekal, Physics Reports , 281 (2001).[6] V. Lucarini, Journal of Statistical Physics , 1698 (2018).[7] Z. Chen, J. de Gier, I. Hiki, and T. Sasamoto, Phys. Rev. Lett. , 240601 (2018).[8] I. Brevik, P. Parashar, and K. V. Shajesh, Phys. Rev. A ,032509 (2018).[9] F. Xu and J. Wang, Phys. Rev. B , 245430 (2014).[10] E. Lenzi, R. Mendes, and A. Rajagopal, Physica A: StatisticalMechanics and its Applications , 503 (2000).[11] F. Cornu and B. Jancovici, J. Chem. Phys. , 2444 (1989).[12] A. Ferrero and G. T´ellez, Journal of Statistical Physics , 759(2007).[13] A. Ferrero and G. T´ellez, Journal of Statistical Mechanics: Theoryand Experiment , P11021 (2014).[14] Y. Nomura, A. S. Darmawan, Y. Yamaji, and M. Imada, Phys.Rev. B , 205152 (2017).[15] D. S. P. Salazar, Phys. Rev. E , 022131 (2017).[16] K. S. Novoselov, A. K. Geim, S. V. Morozov,D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grig-orieva, and A. A. Firsov, Science , 666 (2004),https://science.sciencemag.org/content/306/5696/666.full.pdf.[17] G. Fiori, F. Bonaccorso, G. Iannaccone, T. Palacios, D. Neumaier,A. Seabaugh, S. Banerjee, and L. Colombo, Nature Nanotechnol-ogy , 768 (2014).[18] F. Schwierz, J. Pezoldt, and R. Granzner, Nanoscale , 8261(2015).[19] R. Sterling, H. Rattanasonti, S. Weidt, K. Lake, P. Srinivasan, S. Webster, M. Kraft, and W. Hensinger, Nature communications , 3637 (2014).[20] C. Flindt, N. A. Mortensen, and A.-P. Jauho, Nano Letters , 2515(2005), pMID: 16351206, https://doi.org/10.1021/nl0518472.[21] F. Sugino, JHEP (2004), 10.1088/1126-6708/2004/03/067.[22] R. G. Sarav´ı, F. Schaposnik, and J. Solomin, Nuclear Physics B , 239 (1981).[23] A. Perera and T. Urbic, Journal of Molecular Liquids (2018),10.1016/j.molliq.2018.05.133.[24] M. F. Ra˜nada and M. Santander, Journal of Mathematical Physics , 5026 (1999), https://doi.org/10.1063/1.533014.[25] E. G. Kalnins, J. M. Kress, and W. Miller, Journal of Mathemati-cal Physics , 053509 (2005), https://doi.org/10.1063/1.1897183.[26] J. M. Speight, Phys. Rev. D , 3830 (1997).[27] D. Mitrea and I. Mitrea, Communications inPartial Differential Equations , 304 (2010),https://doi.org/10.1080/03605302.2010.489629.[28] D. G. Truhlar, Journal of Computational Physics , 123 (1972).[29] Z. Jomaa and C. Macaskill, Journal of Computational Physics ,488 (2005).[30] J. L. Steger, Computer Methods in Applied Mechanics and Engi-neering , 175 (1978).[31] T. Ozis, E. Aksan, and A. ¨Ozde¸s, Applied Mathematics and Com-putation , 417 (2003).[32] Y. Lin and C. Xu, Journal of Computational Physics , 1533(2007).[33] J. Izadian, N. Ranjbar, and M. Jalili, World Applied SciencesJournal , 95 (2013).[34] G. Jo and D. Kwak, Numerical Algorithms (2018),10.1007/s11075-018-0544-9.[35] D. Y. Kwak, H. J. Kwon, and S. Lee, Applied Mathematics andComputation , 77 (1999).[36] B. Fornberg, Mathematics of Computation , 699 (1988).[37] G. Forsythe and W. Wasow, Finite Difference Methods for Par-tial Differential Equations: Applied Mathematics Series (Liter-ary Licensing, LLC, 2013). Appendix A: Deduction of the weight function
The relation obeyed by the weight function that makes the Liouville operator self-adjoint can be deduced byperforming a direct substitution of eq. (1) into eq. (5) and using Green’s: R V φ ( ∇ ψ ) d r = R V ψ ( ∇ φ ) d r + H ∂ V h φ ( ~ ∇ ψ ) − ψ ( ~ ∇ φ ) i · n d S and the Divergence: R V a · ( ~ ∇ ψ ) d r = H ∂ V ψ a · n d S − R V ψ ( ~ ∇· a ) d r theorems. After writing it conveniently,the result of this operation is ψ ( r ) = Z r ′ G ( r ′ , r ) h w ( r ′ , r ) (cid:2) ∇ { r ′ } ψ ( r ′ ) (cid:3) + (cid:2) ~ ∇ { r ′ } w ( r , r ′ ) (cid:3) · (cid:2) ~ ∇ { r ′ } ψ ( r ′ ) (cid:3) + w ( r ′ , r ) g ( r ′ ) ψ ( r ′ ) i d r ′ + Z r ′ n(cid:2) ~ ∇ { r ′ } w ( r ′ , r ) − w ( r ′ , r ) ~f ( r ′ ) (cid:3) · ~ ∇ { r ′ } ψ ( r ′ ) + ~ ∇ { r ′ } · (cid:2) ~ ∇ { r ′ } w ( r ′ , r ) − w ( r ′ , r ) ~f ( r ′ ) (cid:3) ψ ( r ′ ) o d r ′ + I ∂ r ′ h w ( r ′ , r ) (cid:2) ψ ( r ′ ) ~ ∇ { r ′ } G ( r ′ , r ) − G ( r ′ , r ) ~ ∇ { r ′ } ψ ( r ′ ) (cid:3) − G ( r ′ , r ) ψ ( r ′ ) (cid:2) ~ ∇ { r ′ } w ( r ′ , r ) − w ( r ′ , r ) ~f ( r ′ ) (cid:3)i · n d S ′ . Notice that under the choice ~ ∇ { r ′ } w ( r ′ , r ) − w ( r ′ , r ) ~f ( r ′ ) = 0, last equation transforms into eq. (7). Appendix B: Finite elements method, matrix elements
The elements of matrices A ( n ) introduced in section III A depend on the required level of accuracy and the centralsite η that we choose; a general algorithm to deduce such elements is shown in [36]. In the simplest case, as we choose η = 0 as the central point in conjunction with the two closets neighbors—three-point stencil approximation—we havethe following relations [36, 37] A (1) = · · · − · · · − · · · ... ... ... . . . ... ... ... · · · · · · − · · · − N +1 × N +1 , A (2) = − · · · − · · · − · · · ... ... ... . . . ... ... ... · · · − · · · − · · · − N +1 × N +1 . (B1)Notice how the matrices A ( n ) must be truncated at the boundaries; this is a natural consequence of the FEM, comingfrom the boundary conditions. Appendix C: Derivatives at the boundaries for DBC
If the Green function is used to find a nonhomgeneous function with DBC, the derivatives of the Green function atthe borders are needed—see eq. (7) and eq. (13). Combining eq. (35) with the boundary conditions stated in eqs. (31)to (34), we deduce the following two relations ( j = 0 and j = N refer to the two possible boundaries) G ′ (0 ,N ) ,kθ,θ ′ = 12 π X λ,µ e iλθ e − iµθ ′ B (0 ,N ) ,kλ,µ , where B (0 ,N ) ,kλ,µ = ∓ P ,N X ν (cid:2) S (0 ,N ) (cid:3) − λ,ν A (1 ,N − ,kν,µ . (C1)The matrix elements associated with S (0 ,N ) are S (0 ,N ) λ,λ = P ,N ∓ h Q ,N and S (0 ,N ) λ,λ − µ = ∓ h Q ,Nµ . Since ψ is known atthe boundaries, the derivatives G ′ , θ,θ ′ and G ′ N,Nθ,θ ′ are irrelevant; additionally, a disk only requires the calculation of G ′ N,kθ,θ ′ . The symmetric elements G ′ j (0 ,N ) θ,θ ′ can be found similarly, in terms of the transpose elements B j, (0 ,N ) λ,µ , whichare defined according to last expression by performing the index change and transposition A (1 ,N − ,kν,µ → A j, (1 ,N − ν,µ .In one dimension the derivatives are G ′ (0 ,N ) k = ∓ P ,N A (1 ,N − ,k P ,N ∓ h Q ,N .8 Appendix D: Symmetry properties of some matrix elements
Expanding eq. (38) to eliminate the negative modes, we can write the Green function as G j,kθ,θ ′ = − h π (cid:20) A j,k , + X λ ≥ n(cid:2) A j,k − λ, + A j,kλ, (cid:3) cos( λθ ) + i (cid:2) − A j,k − λ, + A j,kλ, (cid:3) sin( λθ ) o + X µ ≥ n(cid:2) A j,k , − µ + A j,k ,µ (cid:3) cos( µθ ′ )+ i (cid:2) A j,k , − µ − A j,k ,µ (cid:3) sin( µθ ) o + X λ, µ ≥ n(cid:2) A j,k − λ, − µ + A j,k − λ,µ + A j,kλ, − µ + A j,kλ,µ (cid:3) cos( λθ ) cos( µθ ′ )+ i (cid:2) A j,kλ, − µ + A j,kλ,µ − A j,k − λ, − µ − A j,k − λ,µ (cid:3) sin( λθ ) cos( µθ ′ ) + i (cid:2) A j,k − λ, − µ − A j,k − λ,µ + A j,kλ, − µ − A j,kλ,µ (cid:3) cos( λθ ) sin( µθ ′ )+ (cid:2) A j,k − λ, − µ − A j,k − λ,µ − A j,kλ, − µ + A j,kλ,µ (cid:3) sin( λθ ) sin( µθ ′ ) o(cid:21) . (D1)Since the Green function must be real for real Liouville operators, we demand that the imaginary contributions oflast expression must vanish. Hence, we have the restrictions stated in eq. (39). Appendix E: Expansion of the Green function as sines and cosines
This expansion allows us to write the Green function as a sum of real elements, explicitly showing that the Greenfunction is real. Using the properties stated in eq. (39) into eq. (D1), we find that G jkθθ ′ = − h π Re (cid:0) A j,k , (cid:1) − hπ X λ ≥ h Re (cid:0) A j,kλ, (cid:1) cos( λθ ) + Re (cid:0) A j,k ,λ (cid:1) cos( λθ ′ ) − Im (cid:0) A j,kλ, (cid:1) sin( λθ ) + Im (cid:0) A j,k ,λ (cid:1) sin( λθ ′ ) i − hπ X λ, µ ≥ h Re (cid:0) A j,kλ,µ (cid:1) cos( λθ − µθ ′ ) + Re (cid:0) A j,kλ, − µ (cid:1) cos( λθ + µθ ′ ) i + hπ X λ, µ ≥ h Im (cid:0) A j,kλ,µ (cid:1) sin( λθ − µθ ′ ) + Im (cid:0) A j,kλ, − µ (cid:1) sin( λθ + µθ ′ ) i . (E1)In the presence of angular symmetry, A jkλµ = A jkλµ δ λµ , so last equation reduces to G jkθθ ′ = − h π Re (cid:0) A j,k , (cid:1) − hπ X λ ≥ h Re (cid:0) A j,kλ,λ (cid:1) cos[ λ ( θ − θ ′ )] − Im (cid:0) A j,kλ,λ (cid:1) sin[ λ ( θ − θ ′ )] i . (E2) Appendix F: Computation of ψ ( r ) as an exponential expansion In this section we will derive expressions for eq. (7) and eq. (13) for both DBC and NBC. Although both approachesmust lead to the same results, it is worthwhile to show how both relations can be found through the formalism wehave described.For convenience, we will split ψ ( r, θ ) into a volume ( V ) and surface ( S ) contribution —the volume contributionis the term containing the integral over r ′ in eq. (7) and eq. (13); the surface contribution is the one containing theintegral over the closed surface ∂ r ′ in the same equations. For DBC and NBC, ψ can be written in discrete coordinatesas ( ψ jθ ) DBC = ( ψ jθ ) V + ( ψ jθ ) DBC S , (F1a)( ψ jθ ) NBC = ( ψ jθ ) V + ( ψ jθ ) NBC S . (F1b)There are many ways to evaluate numerically an integral. We will use one of the simplest, however, very efficient, waysto do so, the so called trapezoid rule. Due to the discretization we have used, this rule will be applied to evaluate thethe radial integrals, appearing in the volume contributions; the integrals over angular coordinates will be evaluateddirectly using the Fourier expansions of the functions involved. Using the weight function φ ( r ′ ) → φ kθ ′ = P λ φ kλ e iλθ ′ , the weight function: w ( r ′ , r ) → w kθ ′ w ( r j ,θ ) = w ( r j ,θ ) P λ w kλ e iλθ ′ —and something similar forthe boundary conditions ψ (0 ,N ) θ ′ and ψ ′ (0 ,N ) θ ′ . Now, we will define the function ξ w ( k, M k,j , η k , θ ) = 12 π Z π dθ ′ X λ, µ e iλθ ′ M k,jλ,µ e − iµθ X ν η kν e iνθ ′ X ρ w kρ e iρθ ′ = X λ, µ, ν e − iµθ M k,jλ,µ η kν w k − λ − ν . (F2)This definition will be used to define the volume - and surface -terms.Since the volume -term can be written as R r N r r ′ dr ′ R π w ( r ′ , r ) G ( r ′ , r ) φ ( r ′ ) dθ ′ ( r = R int , r N = R ext ), we can saythat ( ψ jθ ) V = − h w ( r j , θ ) h N − X k =1 r k ξ w ( k, A k,j , φ k , θ ) + 12 r ξ w (0 , A ,j , φ , θ ) + 12 r N ξ w ( N, A N,j , φ , θ ) i . (F3)The surface –term that arises in DBC can be expanded as r ′ R π w ( r ′ , r ) ψ ( r ′ ) ∂ r ′ G ( r ′ , r ) dθ ′ (cid:12)(cid:12) r N r . Similarly as shownabove, in discrete coordinates it is given by( ψ jθ ) DBC S = 1 w ( r j , θ ) h r N ξ w ( N, B N,j , ψ N , θ ) − r ξ w (0 , B ,j , ψ , θ ) i . (F4)Finally, the surface –term − H ∂ r ′ w ( r ′ , r ) G ( r ′ , r ) ~ ∇ { r ′ } ψ ( r ′ ) · n dS ′ that appears in NBC is now expanded as − r ′ R π w ( r ′ , r ) G ( r ′ , r ) ∂ r ′ ψ ( r ′ ) dθ ′ (cid:12)(cid:12) r N r , it now becomes( ψ jθ ) NBC S = hw ( r j , θ ) h r N ξ w ( N, A N,j , ψ N , θ ) − r ξ w (0 , A ,j , ψ , θ ) i . (F5)Remarks: the matrix elements of matrices A and B (0 ,N ) are given by eq. (38) and eq. (C1), respectively—the indicesassociated to the position in the blocks have been omitted by convenience. The function ψ jθ is defined in the interval1 ≤ j ≤ N −
1; in DBC the terms ψ θ and ψ Nθ are given, in NBC the function at the borders is not accurate enoughdue to the discontinuity of the Green function at the borders. Using no weight function
When we adapt the convention stated in eq. (13), eqs. (F3)–(F5) are slightly modified. We now define the function ξ as ξ ( k, M j,k , η k , θ ) = Z π dθ ′ π X λ, µ e iλθ M j,kλ,µ e − iµθ ′ X ν η kν e iνθ ′ = X λ, µ e iλθ M j,kλ,µ η kµ . (F6)We can now conclude that( ψ jθ ) V = − h N − X k =1 h r k ξ ( k, A j,k , φ k , θ ) + 12 r ξ (0 , A j, , φ , θ ) + 12 r N ξ ( N, A j,N , φ N , θ ) i (F7)( ψ jθ ) DBC S = r N ξ ( N, B j,N , ψ N , θ ) − r ξ (0 , B j, , ψ , θ ) (F8)( ψ jθ ) NBC S = hr N (cid:2) ξ ( N, A j,N , f Nr , θ ) + ξ ( N, A j,N , ψ ′ N , θ ) (cid:3) − hr (cid:2) ξ (0 , A j, , f r , θ ) + ξ (0 , A j, , ψ ′ , θ ) (cid:3) . (F9) Appendix G: Computation of the inhomogeneous function as expansion of trigonometric functions
It is now useful to expand the relations shown in previous section as trigonometric functions. Although theexpressions found are much longer, this allows us to use the symmetry properties, eq. (39), to get rid of irrelevantterms and explicitly express ψ ( r ) as a real function. Besides, the exponential expansion defined in appendix F mightintroduce some spurious imaginary contributions, which may arise by as a consequence of the truncating process ofmatrix U —the complex conjugate counterparts of some modes may be discarded in this process. Taking advantage ofthe definitions used in appendix F, eq. (F3) to eq. (F5) are still valid when we adopt the convention stated in eq. (7);0similarly, when the convention eq. (13) is adopted, eq. (F7) to eq. (F9) are also valid. Now, we only need to expand ξ w and ξ eliminating the negative complex modes to express them as sum of real modes. By doing so, eq. (F2) becomes ξ w ( k, M k,j , η k , θ ) = Re (cid:0) M k,j , (cid:1) Re (cid:0) η k (cid:1) Re (cid:0) w k (cid:1) + 2 X λ ≥ Re (cid:0) η k (cid:1)h Re (cid:0) M k,j , (cid:1) X kλ + Y k,jλ + Re (cid:0) w k (cid:1)(cid:2) Re (cid:0) M k,j ,λ (cid:1) cos( λθ ) + Im (cid:0) M k,j ,λ (cid:1) sin( λθ ) (cid:3)i + 4 X λ, µ ≥ X kµ (cid:2) Re (cid:0) M k,j ,λ (cid:1) cos( λθ ) + M I k,j ,λ sin( λθ ) (cid:3) + 2 X λ, µ ≥ Re (cid:0) η k (cid:1)h(cid:2) Re (cid:0) w kλ (cid:1) M (1) k,j (+) λ,µ + Im (cid:0) w kλ (cid:1) M (2) k,j (+) λ,µ (cid:3) cos( µθ ) − (cid:2) Im (cid:0) w kλ (cid:1) M (1) k,j ( − ) λ,µ − Re (cid:0) w kλ (cid:1) M (2) k,j ( − ) λ,µ (cid:3) sin( µθ ) i + 2 X λ,µ ≥ h Re (cid:0) η kµ (cid:1)(cid:2) Re (cid:0) M k,jλ, (cid:1) w (1) k +( λ,µ ) + Im (cid:0) M k,jλ, (cid:1) w (2) k +( λ,µ ) (cid:3) + Im (cid:0) η kµ (cid:1)(cid:2) Re (cid:0) M k,jλ, (cid:1) w (2) k − ( λ,µ ) − Im (cid:0) M k,jλ, (cid:1) w (1) k − ( λ,µ ) (cid:3)i + 2 X λ,µ,ν ≥ h M (1) k,j (+) λ,µ (cid:2) Z (1) k ( λ,ν ) + Z (2) k ( λ,ν ) (cid:3) + M (2) k,j (+) λ,µ (cid:2) Z (3) k ( λ,ν ) + Z (4) k ( λ,ν ) (cid:3)i cos( µθ )+ 2 X λ,µ,ν ≥ h M (2) k,j ( − ) λ,µ (cid:2) Z (1) k ( λ,ν ) + Z (2) k ( λ,ν ) (cid:3) − M (1) k,j ( − ) λ,µ (cid:2) Z (3) k ( λ,ν ) − Z (4) k ( λ,ν ) (cid:3)i sin( µθ ) , (G1)where we used the definitions X kλ = Re (cid:0) η kλ (cid:1) Re (cid:0) w kλ (cid:1) + Im (cid:0) η kλ (cid:1) Im (cid:0) w kλ (cid:1) , Y k,jλ = Re (cid:0) M k,jλ, (cid:1) Re (cid:0) w kλ (cid:1) + Im (cid:0) M k,jλ, (cid:1) Im (cid:0) w kλ (cid:1) ; M (1) k,j ( ± ) λ,µ = Re (cid:0) M k,jλ,µ (cid:1) ± Re (cid:0) M k,jλ, − µ (cid:1) , M (2) k,j ( ± ) λ,µ = Re (cid:0) M k,jλ,µ (cid:1) ± Re (cid:0) M k,jλ, − µ (cid:1) ; w (1) k ± ( λ,ν ) = Re (cid:0) w kλ + ν (cid:1) ± Re( w kλ − ν (cid:1) , w (2) k ± ( λ,ν ) = Im (cid:0) w kλ + ν (cid:1) ± Im( w kλ − ν (cid:1) ; Z (1) k ( λ,ν ) = Re (cid:0) η kν (cid:1) w (1) k +( λ,ν ) , Z (2) k ( λ,ν ) = Im (cid:0) η kν (cid:1) w (2) k − ( λ,ν ) , Z (3) k ( λ,ν ) = Re (cid:0) η kν (cid:1) w (2) k +( λ,ν ) , Z (4) k ( λ,ν ) = Im (cid:0) η kν (cid:1) w (1) k − ( λ,ν ) . (G2)Similarly, eq. (F6) is now written as ξ ( k, M j,k , η k , θ ) = Re (cid:0) M j,k , (cid:1) Re (cid:0) η k (cid:1) + 2 X λ ≥ h R j,kλ + Re (cid:0) η k (cid:1)(cid:2) Re (cid:0) M j,kλ, (cid:1) cos( λθ ) − Im (cid:0) M j,kλ, (cid:1) sin( λθ ) (cid:3)i + 2 X λ, µ ≥ (cid:2) M (1) j,k (+) λ,µ Re (cid:0) η kµ (cid:1) − M (2) j,k ( − ) λ,µ Im (cid:0) η kµ (cid:1)(cid:3) cos( λθ ) − X λ, µ ≥ (cid:2) M (2) j,k (+) λ,µ Re (cid:0) η kµ (cid:1) + M (1) k,j ( − ) λ,µ Im (cid:0) η kµ (cid:1)(cid:3) sin( λθ ) , (G3)with R j,kλ = Re (cid:0) M j,k ,λ (cid:1) Re (cid:0) η kλ (cid:1) − Im (cid:0) M j,k ,λ (cid:1) Im (cid:0) η I kλ (cid:1) . The one dimensional case
The function ψ , satisfying the equation L x ψ ( x ) = φ ( x ), where L x is defined according to eq. (40), can be found bymeans of the relations below. For DBC with either weight function or not ψ j = − h h N − X k =1 w k,j A k,j φ k + 12 w ,j A ,j φ + 12 w N,j A N,j φ N i + P N w N,j ψ N A N − ,j h Q N + P w ,j ψ A ,j − h Q ; (G4a)= − h h N − X k =1 A j,k φ k + 12 A j, φ + 12 A j,N φ N i + q N ψ N A j,N − h Q N + P ψ A j, − h Q . (G4b)For NBC ψ j = − h h N − X k =1 w k,j A k,j φ k + 12 w ,j A ,j φ + 12 w N,j A Nj φ N i + hP k w k,j ψ ′ k A k,j (cid:12)(cid:12)(cid:12) Nk =0 ; (G5a)= − h h N − X k =1 A j,k φ k + 12 A j, φ + 12 A j,N φ N i + hP k A j,k (cid:2) ψ ′ k + ψ k f k (cid:3)(cid:12)(cid:12)(cid:12) Nk =0 ..