Approximate Solution of Two Dimensional Disc-like Systems by One Dimensional Reduction: An Approach through the Green Function Formalism Using the Finite Elements Method

Abstract

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 of an effective one dimensional system. After elaborating on the formalism, we propose to secure an exact solution via a Fourier expansion of the Green function, which entails solving an infinitely countable system of differential equations for the Green–Fourier modes that in the simplest case yields the source-free Green distribution. We present results on non separable systems—or such whose solution cannot be obtained by the usual variable separation technique—on both annulus and disc geometries, and show how the resulting one dimensional Fourier modes potentially generate a near-exact solution. Numerical solutions will be obtained via finite differentiation using Finite Difference Method (FDM) or Finite Element Method (FEM) with the three-point stencil approximation to derivatives. Comparing to known exact solutions, our results achieve an estimated numerical relative error below ${10}^{-6}$. Solutions show the well-known presence of peaks when $\mathbf{r}={\mathbf{r}}^{\prime }$ and a smooth behavior otherwise, for differential equations involving well-behaved functions. We also verified how the Green functions are symmetric under the presence of a “weight function”, which is guaranteed to exist in the presence of a curl-free vector field. Solutions of non-homogeneous differential equations are also shown using the Green formalism and showing consistent results.

1. Introduction

In the present work, we elaborate on the FEM for solving complex two-dimensional partial differential equations (DE) using a Green function construction. Green’s method has been employed extensively in physics for solving Laplace’s equation and associates in a cornucopia of areas, such as Quantum and Statistical Mechanics. In quantum mechanics, for example, the method of non-equilibrium Green’s functions (NEGF) has been used to study the Brownian motion of a quantum oscillator [1], quantum thermal transport [2,3], electron transport through single and two levels of interacting quantum dot [4], derive quantum kinetic equations [5,6], study hadronic physics [7], among others. In statistical mechanics, some of the applications of the Green functions include the predictions of some observables [8], help to describe 1D hydrodynamic models [9], finding electrical properties of some physical systems [10,11], study non-extensive statistical mechanics with new normalized q-expectation values [12], as well as time- and ensemble-average statistical mechanics of the Gaussian network model [13], and so much more. Even the Green functions are used in quantum field theory to describe the propagators of quantum fields in the perturbative regime.

Not only are Green functions useful to solve systems described by inhomogeneous differential equations, but they can also be used to describe thermodynamic properties. For instance, the density and correlations of particles immersed in two-dimensional two component plasmas at certain temperatures can be described by sets of Green functions [14,15,16]. The calculation of the energy gap in the one-dimensional Hubbard model can also be approached by means of the Green function formalism [17]. From a mathematical and more fundamental perspective, the Green functions and domains with uniformly rectifiable boundaries of all dimensions have been approximated [18]; an analysis of extra critical points of Green functions on Flat Tori has also been performed, in which their minimality has been emphazised [19].

In order to study how an inhomogenous partial DE can be solved by the method we propose, we start defining a differential operator

(1)$\begin{array}{c}\hfill {\widehat{\mathcal{L}}}_{\left\{\mathbf{r}\right\}}\square =({\overrightarrow{\nabla }}_{\left\{\mathbf{r}\right\}}+\overrightarrow{f}\left(\mathbf{r}\right))·({\overrightarrow{\nabla }}_{\left\{\mathbf{r}\right\}}\square )+g\left(\mathbf{r}\right)\square ,\end{array}$

acting on a scalar field, denoted by □, in ${\mathbb{R}}^d$, with d the dimension of the system—i.e., $\mathbf{r}\in {\mathbb{R}}^d$. This operator is known in other contexts as the Liouville operator; via this definition, we often describe the evolution of a relevant quantity $\psi$ by means of the equation ${\partial }_t\psi (\mathbf{r},t)-\widehat{\mathcal{L}}\psi (\mathbf{r},t)=0$ as it is the case of the wave function in quantum mechanics. For example, in the diffusion phenomenon the functions take the form $\overrightarrow{f}(\mathbf{r},t)=\overrightarrow{\nabla }D(\mathbf{r},t)$ and $g\left(\mathbf{r}\right)=0$, and for the Helmholtz equation $\overrightarrow{f}(\mathbf{r},t)=0$ and $g\left(\mathbf{r}\right)=m^2$, with m a constant.

Finding solutions to the latter has motivated the development of numerical methods that grow in number and complexity. For instance, using Restricted Boltzmann Machines we can engineer an artificial neural network that is able to accurately sample the probability distribution for quantum statistical systems [20,21]. However, some effort can be made from a mathematical point of view prior to implementing a full scale numerical calculation.

Green’s function—or more precisely, distribution—is perhaps the most interesting artifact of a huge bag of tricks that we have when facing differential equations. Its power relies on the possibility of inverting the differential operator $\widehat{\mathcal{L}}$ to solve the inhomogeneous equation

(2)$\begin{array}{c}\hfill \widehat{\mathcal{L}}\psi \left(\mathbf{r}\right)=\varphi \left(\mathbf{r}\right),\end{array}$

with $\psi \left(\mathbf{r}\right)$ and $\varphi \left(\mathbf{r}\right)$ two scalar functions. Hinting that its existence, the Green Distribution, is conditioned by some properties of $\widehat{\mathcal{L}}$.

A disadvantage of the Green methodology is the duplication of degrees of freedom, encouraging researchers to find $\psi \left(\mathbf{r}\right)$ directly. Our aim is not to develop a generalized theory for an arbitrary problem and number of dimensions. Despite this, we can look into the consequences of breaking down one dimension by focusing on the simple two-dimensional case.

Two-dimensional systems are of great interest in statistical mechanics [14,15,16], material sciences [22,23,24], quantum computing [25,26], high-energy physics [27,28], ionic fluids [29], theoretical mathematics [30,31], and many others.

The outline of the paper is as follows. We first remind some relevant known results for the Green’s function construction in Section 2 prior to presenting the strategy to move from 2D to 1D in Section 2.4. We lay out a clever geometric interpretation of the result in Section 2.1 followed by a connection to a relevant theory for Hilbert Space functions in Section 2.3. Consequently, we next discuss its implications towards finding the Green function using FDM in Section 3. We then present some mathematical results that include the solution of certain well known systems for testing purposes and the implementation of the method in a non-separable 2D system in Section 4; a discussion of such results and an explanation of how the algorithm can be adapted to solve the heat diffusion problem in thermal equilibrium is stated in Section 5. Finally, we wrap up the conclusions in Section 6. Intermediate calculations and numerical details are left for further inspection in appendices.

2. Framework and Methodology

We start studying the Green function formalism by postulating the convolution identity from the Dirac distribution,

(3)$\begin{array}{c}\hfill \psi \left(\mathbf{r}\right)={\int }_{{\mathbf{r}}^{\prime }}\psi \left({\mathbf{r}}^{\prime }\right)\delta ({\mathbf{r}}^{\prime }-\mathbf{r})w({\mathbf{r}}^{\prime },\mathbf{r})\mathrm{d}{\mathbf{r}}^{\prime },\end{array}$

with $w({\mathbf{r}}^{\prime },\mathbf{r})$ a weight function properly defined by two conditions; the first of which $w(\mathbf{r},\mathbf{r})=1$. Now by defining $G({\mathbf{r}}^{\prime },\mathbf{r})$ as,

(4)$\begin{array}{c}\hfill {\widehat{\mathcal{L}}}_{\left\{{\mathbf{r}}^{\prime }\right\}}G({\mathbf{r}}^{\prime },\mathbf{r})=\delta ({\mathbf{r}}^{\prime }-\mathbf{r}),\end{array}$

with $\delta ({\mathbf{r}}^{\prime }-\mathbf{r})=\delta (\mathbf{r}-{\mathbf{r}}^{\prime })$ the ${\mathbb{R}}^d$ Dirac delta distribution, then,

(5)$\begin{array}{c}\hfill \psi \left(\mathbf{r}\right)={\int }_{{\mathbf{r}}^{\prime }}\psi \left({\mathbf{r}}^{\prime }\right)\left[{\widehat{\mathcal{L}}}_{\left\{{\mathbf{r}}^{\prime }\right\}}G({\mathbf{r}}^{\prime },\mathbf{r})\right]w({\mathbf{r}}^{\prime },\mathbf{r})\mathrm{d}{\mathbf{r}}^{\prime }.\end{array}$

The second condition over $w(\mathbf{r},{\mathbf{r}}^{\prime })$ will be determined in such a way that $\widehat{\mathcal{L}}$ is self-adjoint (Hermitian), or equivalently

(6)$\begin{array}{cc}\hfill {\int }_{{\mathbf{r}}^{\prime }}\psi \left({\mathbf{r}}^{\prime }\right)& \hfill \left[{\widehat{\mathcal{L}}}_{\left\{{\mathbf{r}}^{\prime }\right\}}G({\mathbf{r}}^{\prime },\mathbf{r})\right]w({\mathbf{r}}^{\prime },\mathbf{r})\mathrm{d}{\mathbf{r}}^{\prime }={\int }_{{\mathbf{r}}^{\prime }}\left[{\widehat{\mathcal{L}}}_{\left\{{\mathbf{r}}^{\prime }\right\}}\psi \left({\mathbf{r}}^{\prime }\right)\right]G({\mathbf{r}}^{\prime },\mathbf{r})w({\mathbf{r}}^{\prime },\mathbf{r})\mathrm{d}{\mathbf{r}}^{\prime }+\mathrm{b}.\mathrm{c}.,\hfill \end{array}$

with added Dirichlet or Neumann boundary conditions (b.c.). Direct substitution into Equation (5), using Equation (2), yields

(7)$\begin{array}{c}\hfill \psi \left(\mathbf{r}\right)={\int }_{{\mathbf{r}}^{\prime }}G({\mathbf{r}}^{\prime },\mathbf{r})\varphi \left({\mathbf{r}}^{\prime }\right)w({\mathbf{r}}^{\prime },\mathbf{r})\mathrm{d}{\mathbf{r}}^{\prime }+\mathrm{b}.\mathrm{c}..\end{array}$

2.1. On the Nature of $w(\mathbf{r},{\mathbf{r}}^{\prime })$, and ${\widehat{\mathcal{L}}}^{-1}$

This former known result deserves a more delicate look, particularly, on the existence of the weight function, and how previous solution relates with the usual convolution theorem $\psi \left(\mathbf{r}\right)={\int }_{{\mathbf{r}}^{\prime }}G(\mathbf{r},{\mathbf{r}}^{\prime })\varphi \left({\mathbf{r}}^{\prime }\right)\mathrm{d}{\mathbf{r}}^{\prime }+\mathrm{b}.\mathrm{c}.$ As mentioned, an appropriate choice for the weight function ensures that Equation (5) reproduces Equation (7). This is performed by using the Green’s and Divergence theorem in Equation (5) to perform an integration by parts. After simplifications, we realize that by choosing the weight function, such that ${\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}w({\mathbf{r}}^{\prime },\mathbf{r})-w({\mathbf{r}}^{\prime },\mathbf{r})\overrightarrow{f}\left({\mathbf{r}}^{\prime }\right)=0$ (See Appendix A for further details) we ensure that the operator is self-adjoint! This is essential to Green’s method. Hence, if no weight function exists, we might be forced to use other analytical and/or numerical procedures in order to find $\psi$. For that matter, the range of problems that we aim to analyze is narrowed down to the few ones satisfying the aforementioned condition; despite this, a great many of this subset are of special interest for mathematics and physics.

Assuming $w({\mathbf{r}}^{\prime },\mathbf{r})$ exists we are able to incorporate the premise for Equation (7) yielding exactly,

(8)$\begin{array}{cc}\hfill \psi \left(\mathbf{r}\right)=& \hfill {\int }_{{\mathbf{r}}^{\prime }}G({\mathbf{r}}^{\prime },\mathbf{r})\varphi \left({\mathbf{r}}^{\prime }\right)w({\mathbf{r}}^{\prime },\mathbf{r})\mathrm{d}{\mathbf{r}}^{\prime }+\hfill \\ \hfill & \hfill {\oint }_{\partial {\mathbf{r}}^{\prime }}w({\mathbf{r}}^{\prime },\mathbf{r})\left[\psi \left({\mathbf{r}}^{\prime }\right){\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}G({\mathbf{r}}^{\prime },\mathbf{r})-G({\mathbf{r}}^{\prime },\mathbf{r}){\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}\psi \left({\mathbf{r}}^{\prime }\right)\right]·\mathbf{n}\mathrm{d}{S}^{\prime }.\hfill \end{array}$

Notice how the second term depends on $\psi \left(\mathbf{r}\right)$’s boundary conditions. By choosing identical conditions and trivial values for the Green distribution function at the boundaries ($G=0$ for Dirichlet or $G^{\prime }=0$ for Neumann) we are capable of solving an infinite number of similar boundary value problems.

The discussion for the existence of the weight function can be answered mathematically. Given the relationship required, the weight is defined as $w({\mathbf{r}}^{\prime },\mathbf{r})\equiv exp[-\gamma ({\mathbf{r}}^{\prime },\mathbf{r})]$, yielding ${\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}\gamma ({\mathbf{r}}^{\prime },\mathbf{r})=-\overrightarrow{f}\left({\mathbf{r}}^{\prime }\right)$. Considering that the curl of the gradient of any scalar function is trivial then $\gamma ({\mathbf{r}}^{\prime },\mathbf{r})$ exists if, and only if, $\overrightarrow{\nabla }\times \overrightarrow{f}=0$, which means that $\overrightarrow{f}$ must be a conservative vector field! Within this view, $\gamma ({\mathbf{r}}^{\prime },\mathbf{r})$ represents the scalar potential associated with a force. Anticipating this last restriction, the solution for $\gamma ({\mathbf{r}}^{\prime },\mathbf{r})$ is independent of a path that simply connects $\mathbf{r}$ to ${\mathbf{r}}^{\prime }$ yielding,

(9)$\begin{array}{c}\hfill w({\mathbf{r}}^{\prime },\mathbf{r})\equiv \frac{e^{-\gamma \left({\mathbf{r}}^{\prime }\right)}}{e^{-\gamma \left(\mathbf{r}\right)}}=\frac{1}{w(\mathbf{r},{\mathbf{r}}^{\prime })},\end{array}$

reflecting on the symmetry of the distribution as it will be shown later. Finally, we are ready to define the inverse operator of $\widehat{\mathcal{L}}$ as

(10)$\begin{array}{c}\hfill {\widehat{\mathcal{L}}}_{\left\{\mathbf{r}\right\}}^{-1}\square ={\int }_{{\mathbf{r}}^{\prime }}G({\mathbf{r}}^{\prime },\mathbf{r})\square \left({\mathbf{r}}^{\prime }\right)w({\mathbf{r}}^{\prime },\mathbf{r})\mathrm{d}{\mathbf{r}}^{\prime },\end{array}$

conditioned by the boundary-value problem, which in turn defines $G({\mathbf{r}}^{\prime },\mathbf{r})$ from Equation (4).

There is one last piece of the puzzle to be resolved and it is related to the symmetry of the Green function distribution. Let us evaluate ${\widehat{\mathcal{L}}}_{\left\{\mathbf{r}\right\}}G({\mathbf{r}}^{\prime },\mathbf{r})$—i.e., the operator acting on the second variable. Direct application of ${\widehat{\mathcal{L}}}_{\left\{\mathbf{r}\right\}}$ on Equation (8), using Equation (2),

(11)$\begin{array}{cc}\hfill \varphi \left(\mathbf{r}\right)=& \hfill {\int }_{{\mathbf{r}}^{\prime }}{\widehat{\mathcal{L}}}_{\left\{\mathbf{r}\right\}}\left\{G({\mathbf{r}}^{\prime },\mathbf{r})w({\mathbf{r}}^{\prime },\mathbf{r})\right\}\varphi \left({\mathbf{r}}^{\prime }\right)\mathrm{d}{\mathbf{r}}^{\prime }+{\widehat{\mathcal{L}}}_{\left\{\mathbf{r}\right\}}\left\{\mathrm{b}.\mathrm{c}.\right\},\hfill \end{array}$

hints how this operator appears to work and leads us to anticipate the convolution of a Dirac distribution. Indeed this is true. To clarify, here we exchanged integral and Liouville operators because they are acting on separate variables, and the weight and Green functions (except at ${\mathbf{r}}^{\prime }=\mathbf{r}$) are differentiable.

This conjecture can be proved from the following statement: two separate problems with different boundary values and identical inhomogeneous differential equation—Equation (2)—share the same Green function distribution and satisfy Equation (11); therefore, by comparing equations for any two cases leads to ${\widehat{\mathcal{L}}}_{\left\{\mathbf{r}\right\}}\left\{\mathrm{b}.\mathrm{c}.\right\}=0$ because we can always choose convenient trivial boundary values (i.e., $\mathrm{b}.\mathrm{c}.=0$) in one case. Consequently,

(12)$\begin{array}{c}\hfill {\widehat{\mathcal{L}}}_{\left\{\mathbf{r}\right\}}\left\{G({\mathbf{r}}^{\prime },\mathbf{r})w({\mathbf{r}}^{\prime },\mathbf{r})\right\}\equiv \delta ({\mathbf{r}}^{\prime }-\mathbf{r}),\end{array}$

an equality that bears meaning in the sense of the distributions. These final result unravels the symmetry of the Green distribution function via the weight function, i.e.,

(13)$\begin{array}{c}\hfill G({\mathbf{r}}^{\prime },\mathbf{r})=G(\mathbf{r},{\mathbf{r}}^{\prime })w(\mathbf{r},{\mathbf{r}}^{\prime }).\end{array}$

An interesting question now arises, and it is related to the possibility of using Equation (13) to drop the weight function out of the equation. This operation, with the addition of the relation ${\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}w(\mathbf{r},{\mathbf{r}}^{\prime })=-w(\mathbf{r},{\mathbf{r}}^{\prime })\overrightarrow{f}\left({\mathbf{r}}^{\prime }\right)$, leads to,

(14)$\begin{array}{cc}\hfill \psi \left(\mathbf{r}\right)=& {\int }_{{\mathbf{r}}^{\prime }}G(\mathbf{r},{\mathbf{r}}^{\prime })\varphi \left({\mathbf{r}}^{\prime }\right)\mathrm{d}{\mathbf{r}}^{\prime }+\hfill \\ \hfill & \hfill {\oint }_{\partial {\mathbf{r}}^{\prime }}\left[\psi \left({\mathbf{r}}^{\prime }\right){\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}G(\mathbf{r},{\mathbf{r}}^{\prime })-G(\mathbf{r},{\mathbf{r}}^{\prime })\left[\psi \left({\mathbf{r}}^{\prime }\right)\overrightarrow{f}\left({\mathbf{r}}^{\prime }\right)+{\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}\psi \left({\mathbf{r}}^{\prime }\right)\right]\right]·\mathbf{n}d{S}^{\prime }.\hfill \end{array}$

Notice that for Neumann boundary conditions (NBC), unlike Dirichlet (DBC), both $\psi \left(\mathbf{r}\right)$ and its derivative—at the boundaries—are necessary. Ergo, Equation (14) is inconvenient for NBC unless either $\psi$ vanishes or $\overrightarrow{f}\left(\mathbf{r}\right)=0$. In such a case, it deems necessary to use the version that incorporates weight function.

Actually, the vector field $\overrightarrow{f}$ does not appear in some of the Liouville operators used in physics. For instance, the Green function associated with the electrostatic field satisfies the relation ${\overrightarrow{\nabla }}^{2}G(\mathbf{r},{\mathbf{r}}^{\prime })=-4\pi \delta (\mathbf{r}-{\mathbf{r}}^{\prime })$—the irrelevant factor of $-4\pi$ appears by convenience. The static regime of the Klein Gordon equation—which also leads to the Yukawa potential—also follows a similar behavior, as its associated Green function in 2D is $K_0\left(\mu r\right)$, satisfying the DE $({\overrightarrow{\nabla }}^2-{\mu }^2)G\left(\mathbf{r}\right)=-2\pi \delta \left(\mathbf{r}\right)$ [32]. Clearly, $\overrightarrow{f}$ is absent in both systems.

Yet, the DEs describing the behavior of other physical systems, such as the driven damped harmonic oscillator, (the one-dimensional driven damped harmonic oscillator is modeled by the DE $m\ddot{x}+b\dot{x}+kx=F\left(t\right)$. The damping constant b plays the role of $\overrightarrow{f}$ in this one dimensional system. Although the time t is the relevant variable describing this system (instead of the position x), the one dimensional formalism we describe is analog to this model.) the diffusion equation at thermal equilibrium with an anisotropic diffusion coefficient, and the electrostatic potential in the presence of anisotropic media, include the existence of a vector field $\overrightarrow{f}$—see Section 4 for more details about the first system.

Surprisingly, any dependence on the weight function in Equation (14) has vanished. As previously stated, the weight function can only be defined when $\overrightarrow{f}$ is a conservative field. Then, an important question now arises: is Equation (14) still valid for non-conservative vector fields? This in a fundamental question that can be addressed in a future work. Since the main purpose is to present a compact and rigorous algorithm to solve the Green function in 2D space, we will restrict our analysis to only the supported cases.

2.2. Boundary Conditions

As previously stated, the Green function conveniently inherits identical types of conditions as the target function $\psi$ at the boundaries. These can be summarized as,

(15)$\begin{array}{cc}\hfill \mathrm{DBC}:\to & \hfill G(\mathbf{r},{\mathbf{r}}^{\prime }){|}_{\mathbf{r}\mathrm{at}{R}_{\mathrm{ext}/\mathrm{int}}}=0,\hfill \end{array}$

(16)$\begin{array}{cc}\hfill \mathrm{NBC}:\to & \hfill {\partial }_{\mathbf{r}}G(\mathbf{r},{\mathbf{r}}^{\prime }){|}_{\mathbf{r}\mathrm{at}{R}_{\mathrm{ext}/\mathrm{int}}}=0.\hfill \end{array}$

However, there are two hidden additional conditions that must be satisfied enforced by the presence of Dirac’s distribution. The rationale behind this is that without them $G=0$ will be a solution to the Green function for the simple boundary value problem. Although this is directly visible for Dirichlet, notice that it also applies for Neumann’s case. The added restrictions appear at the artificial boundary $\mathbf{r}={\mathbf{r}}^{\prime }$ implying continuity of G and discontinuity of the local derivative. Both are essential to secure a non-zero solution. Continuity is often regarded considering that the Green distribution is still a function and its derivatives up to second order exist in the classical sense of the DE everywhere except at $\mathbf{r}={\mathbf{r}}^{\prime }$. Though there is a stronger argument that stems from the fact that the annulus and the disc are Lipschitz domains [33], in DE it always results convenient to decide what do we take as an acceptable solution to any problem, which is our particular case here.

Turning to the plane $\mathbf{r}={\mathbf{r}}^{\prime }$, conditions are derived directly from Equation (4) by integrating over $\mathbf{r}$ inside the volume delimited by the surface $S_{\delta }$ enclosing ${\mathbf{r}}^{\prime }$ such that it is contained inside a vicinity—${\mathcal{V}}_{\delta }$—of ${\mathbf{r}}^{\prime }$ (see right of Figure 1 for an artistic view). Using the divergence theorem, the condition simplifies to,

(17)$\begin{array}{c}\hfill \underset{\delta \to 0}{lim}{\oint }_{S_{\delta }}{\overrightarrow{\nabla }}_{\left\{\mathbf{r}\right\}}G(\mathbf{r},{\mathbf{r}}^{\prime })·\widehat{n}dS=1,\end{array}$

where we have kept the leading contributing term while taking the limit. For the one dimensional case, it yields the relation $G^{\prime }(r_{>}^{\prime },r^{\prime })-G^{\prime }(r_{<}^{\prime },r^{\prime })=1$.

2.3. Connection with the Sturm–Liouville Problem

The weight function, if existent, is able to transform the Liouville operator into a self-adjoint differential operator. Notice that action of $w({\mathbf{r}}^{\prime },\mathbf{r})$ on Equation (1),

(18)$w({\mathbf{r}}^{\prime },\mathbf{r}){\widehat{\mathcal{L}}}_{\left\{{\mathbf{r}}^{\prime }\right\}}\square ={\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}·\left[w({\mathbf{r}}^{\prime },\mathbf{r}){\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}\square \right]+w({\mathbf{r}}^{\prime },\mathbf{r})g\left({\mathbf{r}}^{\prime }\right)\square ,$

yields the otherwise known Sturm–Liouville form for Partial Differential Equations (PDE)’s. Namely, the Sturm–Liouville differential operator reads then as,

(19)$\begin{array}{c}\hfill {\widehat{\mathcal{L}}}_{\left\{{\mathbf{r}}^{\prime }\right\}}^{\mathrm{SL}}\square =w({\mathbf{r}}^{\prime },\mathbf{r}){\widehat{\mathcal{L}}}_{\left\{{\mathbf{r}}^{\prime }\right\}}\square .\end{array}$

Consequently, operating onto the Green distribution function gives equal results for both operators, i.e., ${\widehat{\mathcal{L}}}_{\left\{{\mathbf{r}}^{\prime }\right\}}^{\mathrm{SL}}G({\mathbf{r}}^{\prime },\mathbf{r})=\delta (\mathbf{r}-{\mathbf{r}}^{\prime })$.

This last result connects the Sturm–Liouville problem with null eigenvalues and the Green function distribution problem where the former is the solution to the first strictly when $\mathbf{r}\ne {\mathbf{r}}^{\prime }$ under either Dirithlet or Neumann boundary conditions. Resulting from this, $G({\mathbf{r}}^{\prime },\mathbf{r})$ is continuous everywhere and differentiable at $\mathbf{r}\ne {\mathbf{r}}^{\prime }$; the behavior of its derivative at $\mathbf{r}={\mathbf{r}}^{\prime }$ is dictaminated by the Liouville operator in the domain of the problem and specified by the Dirac Delta distribution.

2.4. Moving from 2D to 1D: The Infinite Coupling

As noted before, let us elaborate on the simplest scenario where a reduction in the dimension of the problem significantly improves our chances of procuring a general solution. Assume we would like to find the two-dimensional Green function in accordance with Equation (4). Taking advantage of the completeness of the Fourier infinite expansion of any periodic function, we propose to solve the two-dimensional DE in polar coordinates; the dimensional reduction occurs due to the periodicity in $\theta$ that does not take place in Cartesian coordinates.

Although the Laplace operator in polar coordinates is known to be separable in the variables r and $\theta$, the introduction of the additional terms in Equation (1), as already mentioned, may lead to a DE that cannot be split conveniently. Equation (4) is then given by,

(20)$\begin{array}{cc}\hfill & \hfill \left[\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}+\frac{1}{r^2}\frac{{\partial }^2}{\partial {\theta }^2}+f_r(r,\theta )\frac{\partial }{\partial r}+f_{\theta }(r,\theta )\frac{1}{r}\frac{\partial }{\partial \theta }+g(r,\theta )\right]G(\mathbf{r},{\mathbf{r}}^{\prime })=\delta (\mathbf{r}-{\mathbf{r}}^{\prime }),\hfill \end{array}$

where convenient periodic conditions that must be satisfied suggests we should expand the Green function distribution in Fourier modes. Tentatively, we can resort to a Fourier expansion of the form ${\sum }_{\lambda }{e}^{i\lambda (\theta -{\theta }^{\prime })}{G}_{\lambda }$, where

$\begin{array}{c}\hfill f(r,\theta )=\sum _{\mu \in \mathbb{Z}}{f}_{\mu }\left(r\right)e^{i\mu \theta };f_{\mu }\left(r\right)=\frac{1}{2\pi }{\int }_0^{2\pi }f(r,\theta )e^{-i\mu \theta }d\theta ,\end{array}$

to match the delta distribution expansion—i.e., $\delta (\mathbf{r}-{\mathbf{r}}^{\prime })=\frac{1}{2\pi r}\delta (r-r^{\prime }){\sum }_{\lambda }{e}^{i\lambda (\theta -{\theta }^{\prime })}$—but since we cannot guarantee that the $G_{\lambda }$ coefficients are ${\theta }^{\prime }$-independent (only under proper angular symmetry conditions) then we will assume $G(\mathbf{r},{\mathbf{r}}^{\prime })$ expands as,

(21)$\begin{array}{c}\hfill G(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{2\pi }\sum _{\lambda \in \mathbb{Z}}{e}^{i\lambda \theta }{G}_{\lambda }(r,r^{\prime },{\theta }^{\prime }).\end{array}$

With that in mind and multiplying Equation (20) by $r^2$ to avoid divergences at $r=0$,

(22)$\begin{array}{cc}\hfill & \hfill \sum _{\lambda \in \mathbb{Z}}{e}^{i\lambda \theta }\left[r^2\frac{d^2}{d{r}^2}+r(1+r{f}_r)\frac{d}{dr}+(-{\lambda }^2+i\lambda r{f}_{\theta }+r^{2}g)\right]G_{\lambda }=\sum _{\lambda \in \mathbb{Z}}{e}^{i(\lambda -{\lambda }^{\prime })\theta }r\delta (r-r^{\prime }).\hfill \end{array}$

Notice how we cannot obtain a solution because there remains a residual dependence of $\theta$ in functions $\overrightarrow{f}$ and g. Despite this, a simplification can be manufactured when they are replaced by their Fourier series form before integrating on $\theta$ over a full period. This step yields our master equation where we deduce that the $G_{\lambda }(r,r^{\prime },{\theta }^{\prime })$ modes satisfy the DE (we have dropped out the dependencies of all functions on r, $\theta$, $r^{\prime }$, and ${\theta }^{\prime }$ facilitating a comprehensible reading),

(23)$\begin{array}{cc}\hfill & r^2{G}_{\lambda }^{″}+r{G}_{\lambda }^{\prime }-{\lambda }^2{G}_{\lambda }+\sum _{\mu \in \mathbb{Z}}{r}^2{f}_{r\mu }{G}_{\lambda -\mu }^{\prime }+\sum _{\mu \in \mathbb{Z}}\left[ir(\lambda -\mu )f_{\theta \mu }+r^2{g}_{\mu }\right]G_{\lambda -\mu }\hfill \\ \hfill & =r\delta (r-r^{\prime })e^{-i\lambda \theta {}^{\prime }}.\hfill \end{array}$

This final result shows we have accomplished to reduce the rank of the effective Green function to solve at the expense of requiring a countable large number of these Green modes. It is remarkable how the dependence on ${\theta }^{\prime }$ is delegated to a quasi-negligible term at the right hand side of the equation. We will develop this argument further in the following sections.

This formulation represents an infinitely coupled system of linear differential equations that unsurprisingly contains the solution to the Green function for the classical source-free wave function; the structure of functions $\overrightarrow{f}$ and g defines the strength of the entanglement of Green’s free wave modes appearing in the rate at which the Fourier coefficients—functions—go to zero with increasing mode frequency. For simplicity, we opt to recall Green’s Fourier modes as $\lambda$-modes, and $\overrightarrow{f}$ and g’s modes as $\mu$-modes suggested by the indexes employed in the equation above.

2.5. More on Boundary Conditions of the $\lambda$-Modes

One last effort must be made to explain how boundary conditions are inherited along the free-wave modes. The key to this understanding depends on the geometry of the problem and the originating expansion from Equation (21); we can identify two cases for disc-like systems: the annulus and the disc. Other geometries will be studied in a future work. For the annulus, either under Dirichlet or Neumann boundary conditions, the function or its derivative must vanish at the boundaries. This can be met if all modes preserve the vanishing values at both inner and outer boundaries—under uniform convergence. In doing so, we guarantee to meet all requisites for the Green function and a solution is obtained. Conversely, preserving boundary conditions for the disc is not trivial because we do not have one but two boundaries (the second at $r\to 0^{+}$). Due to the oscillating behavior of $e^{i\lambda \theta }$ with $\theta$ at $r\to 0^{+}$, all $\lambda \ne 0$ modes must vanish at the origin to ensure continuity of the Green distribution function. This can be enforced examining Equation (23) as $r\to 0^{+}$. Discontinuity due to the source at $r=r^{\prime }$ may be neglected for now to realize that we can, while approaching the origin, consider the behavior of each $r^0$, $r^1$, and $r^2$ terms independently. We draw then conveniently,

(24)$\begin{array}{cc}\hfill r^0\left[-{\lambda }^2{G}_{\lambda }\right]& \hfill \underset{r\to 0}{=}0,\hfill \end{array}$

(25)$\begin{array}{cc}\hfill r^1\left[G_{\lambda }^{\prime }+i\sum _{\mu \in \mathbb{Z}}(\lambda -\mu )f_{\theta \mu }{G}_{\lambda -\mu }\right]& \hfill \underset{r\to 0}{=}0,\hfill \end{array}$

(26)$\begin{array}{cc}\hfill r^2\left[G_{\lambda }^{″}+\sum _{\mu \in \mathbb{Z}}{f}_{r\mu }{G}_{\lambda -\mu }^{\prime }+\sum _{\mu \in \mathbb{Z}}{g}_{\mu }{G}_{\lambda -\mu }\right]& \hfill \underset{r\to 0}{=}0,\hfill \end{array}$

where we can choose, via $r^0$ terms, that $G_{\lambda }=0\forall \lambda \ne 0$ and $G_0\ne 0$. Plugging this sequentially into $r^1$ and $r^2$ terms hints $G_{\lambda }^{\prime }=0$ and $G_{\lambda }^{″}=0$ ($\forall \lambda \ne 0$) assuming that ${lim}_{r\to 0}{f}_{\theta \mu }{G}_{\lambda -\mu }=0$ and ${lim}_{r\to 0}{f}_{r\mu }{G}_{\lambda -\mu }^{\prime }=0\wedge {lim}_{r\to 0}{g}_{\mu }{G}_{\lambda -\mu }=0$ (with the exception of $\lambda =0$ where the $\mu =0$ term remains, thus we will choose $G_{\lambda }^{″}=-g_0{G}_0$), respectively.

This is supported from continuity of $g\left(\mathbf{r}\right)$ everywhere in the disc and from the definition of $\overrightarrow{f}\left(\mathbf{r}\right)$, where limited by the existence of $w({\mathbf{r}}^{\prime },\mathbf{r})$, as the gradient of a scalar function. If such functions, $\gamma \left(\mathbf{r}\right)$, where to be free of pathologies and differentiable everywhere in the disc (including the origin) then ${lim}_{r\to 0}{f}_{r\lambda }=0$ and ${lim}_{r\to 0}{f}_{\theta \lambda }=0$ for $\lambda \ne 0$. It remains to say, that in order to fulfill all above conditions we will require that $f_{\theta 0}$ and $f_{r0}$ are finite as $r\to 0^{+}$. Looking under the hood of these assumptions, note that consequently the 0-mode has a logarithmic divergence when $r^{\prime }=0$, i.e., $G_0\underset{r\to 0}{\propto }-logr$.

In summary, the conditions for the disc at the origin are the following two only for $r^{\prime }>0$ (see Section 3 for numerical details)

(27)$G_{\lambda }(0^{+},r^{\prime },{\theta }^{\prime })=0\forall \lambda \ne 0,$

(28)$G_{\lambda }^{\prime }(0^{+},r^{\prime },{\theta }^{\prime })=0\forall \lambda .$

Exceptions and particularities emerging from the specific form of functions $\overrightarrow{f}$ and g must be taken into account when detailing the boundary conditions and may alter the relationships obtained above.

This relationship has to be completed with the resulting relationship at the artificial boundary $r=r^{\prime }$ obtained when using a complementary surface $S_{ϵ}$ corresponding to an open ring of $ϵ>0$ thickness—see left Figure 1 and Equation (17). This gives,

(29)$\begin{array}{c}\hfill r^{\prime }\underset{ϵ\to 0}{lim}{\int }_0^{2\pi }\left[{\partial }_{r}G{(\mathbf{r},{\mathbf{r}}^{\prime })}_{>}-{\partial }_{r}G{(\mathbf{r},{\mathbf{r}}^{\prime })}_{<}\right]d\theta =1.\end{array}$

which entails the radial averaged contribution. We have eliminated angular contributions by selecting the convenient contour $S_{ϵ}$ suggesting a pathway to extend it to the $\lambda$-modes. Direct substitution of Equation (21) along with a convenient choice of unity—inspired by Equation (23)—gives us ultimately,

(30)$\begin{array}{c}\hfill r^{\prime }\left[G{}_{\lambda }^{\prime }(r_{>}^{\prime },r^{\prime },{\theta }^{\prime })-G{}_{\lambda }^{\prime }(r_{<}^{\prime },r^{\prime },{\theta }^{\prime })\right]=e^{-i\lambda {\theta }^{\prime }},\end{array}$

where primes denote partial derivatives with respect to the first argument at both left (<) and right (>) hand sides of $r^{\prime }$. Substituting this relationship in the differential equation, we obtain a similar relationship for the second derivatives, essential to the numerical method, as follows,

(31)${\left(r^{\prime }\right)}^2[G{}_{\lambda }^{″}(r_{>}^{\prime },r^{\prime },{\theta }^{\prime })-G{}_{\lambda }^{″}(r_{<}^{\prime },r^{\prime },{\theta }^{\prime })]=-e^{-i\lambda {\theta }^{\prime }}\left[1+r^{\prime }{f}_r(r^{\prime },{\theta }^{\prime })\right].$

For the disc, the $r^{\prime }=0$ case must be clarified. In polar coordinates, Dirac’s distribution is best described as absent of angular dependence, which entails that for all $\lambda$-modes except $\lambda =0$ it is exactly zero. Therefore, the boundary at the origin for each $\lambda$-mode is dictated by symmetry except for $\lambda =0$. This last, carries the logarithmic divergence. This means that conditions for non-zero modes are unchanged. For the zero mode and due to symmetry $G_0^{\prime }{{|}_{r\to 0^{+}}=-G_0^{\prime }|}_{r\to 0^{-}}$ and $G_0^{″}{{|}_{r\to 0^{+}}=G_0^{″}|}_{r\to 0^{-}}$; however, due to the divergence a cut-off must be set in place. Such a choice of cut-off will be discussed later.

3. Finite Differences Method, FDM, or FEM on a Regular Grid

The FDM, or uniform mesh FEM, has been used extensively in the literature to find approximate solutions for many physical systems and its stability makes it a suitable candidate to obtain a numerical Green distribution function. Some examples include the one-dimensional Schrödinger equation [34], the Poisson equation for Electrodynamics [35], the Euler equations of inviscid fluid flow [36], solutions to 1D and 2D Burgers’ equation [37], and the time-fractional diffusion equation [38]. From the mathematical perspective, the same method has been implemented to solve elliptic, hyperbolic and parabolic partial DEs on irregular meshes [39], with interfaces [40], or in finding optimal algorithms on non-trivial meshes [41].

Orchestrating an exact solution to Equation (23) is virtually not possible. There are four cases where an analytical approach can be attempted: two cases where either $\overrightarrow{f}$ or g are zero, requiring to find a base of $G_{\lambda }$’s that can decouple the system—hence, a diagonalization—the unique case where the same base applies to both coupling matrices accompanying $G_{\lambda }^{\prime }$ and $G_{\lambda }$, and the trivial free-wave ($\overrightarrow{f}$ and g zero). Excluding the latter, finding this diagonalizing operator for the first three cases will be addressed in a future study.

Therefore, we will compute a numerical solution where we approximate the operator with finite differentiation (the finite difference method—FDM or FEM for a regular grid) and bind expansions to include all relevant Green and function modes up to a calculated cut-off; maximum and minimum modes will be chosen, respectively, for $\lambda$- and $\mu$-modes symmetrically as $\left|\lambda \right|\le L$ and $\left|\mu \right|\le M$ considering that $M\le L$ for reasons that will be clarified afterwards.

Since the Green function is twice differentiable, when $\mathbf{r}\ne {\mathbf{r}}^{\prime }$, its Fourier series converges uniformly and its coefficients decay at least as ${\lambda }^{-2}$, conditioned by equally well-behaved functions $\overrightarrow{f}$ and g. Then, a possible educated choice of L is the minimum integer such that the sum of $1/k^2$ up to L exceeds $\frac{{\pi }^2}{6}p$, with p a percentage of accuracy; for example, to achieve at most $1\%$ of estimation error we require $L>60$.

Numerical details and calculations performed henceforth are presented solely for the 3-point stencil. The strategy for the implementation of more accurate approximations will only be mentioned and briefly discussed; their details will be left for the reader to carry them out. Other minor and major details regarding the procedure will be addressed in future works.

3.1. A Large Matrix Equation

The Finite Differences Method (FDM) or Finite Elements Method (FEM) with uniform grid is a simple approach to computing derivatives of functions at a point by using Taylor expansions on a discretized mesh. (The choice of whether dissecting uniformly or non-uniformly is highly dependable on the problem. For example, if we were interested in fracture dynamics we would prefer a non-uniform grid to model complex material topologies.) In doing so, a derivative will rely on knowledge of the values of the function in neighboring sites. Such is the art of computing derivatives. The number of neighboring sites to be taken into consideration determines the degree of which the function approaches to the point value. For instance, in the so called three-point stencil (the site in question and its two adjacent neighbors), the first and second derivatives are accurate up to order square of the mesh size.

Going back to our problem in Equation (23), we turn to a simply redefined one dimensional DE for a sketch of the forthcoming operations. The left hand side reads rewritten as,

$\begin{array}{c}\hfill P\left(r\right)\frac{d^2}{d{r}^2}{G}_{\lambda }(r,r^{\prime },{\theta }^{\prime })+\sum _{\mu \in \mathbb{Z}}{Q}_{\mu }\left(r\right)\frac{d}{dr}{G}_{\lambda -\mu }(r,r^{\prime },{\theta }^{\prime })+\sum _{\mu \in \mathbb{Z}}{R}_{\lambda ,\mu }\left(r\right)G_{\lambda -\mu }(r,r^{\prime },{\theta }^{\prime }).\end{array}$

In transforming the continuous variables r and $r^{\prime }$ into a discrete equally-spaced mesh of size h, we will adopt matrix notation for variables and functions; ergo, for N partitions defining $N+1$ points $h=(R_{\mathrm{ext}}-R_{\mathrm{int}})/N$ in a disc-like geometry. For clarity, we summarize notation changes in Table 1. This procedure applied over the aforementioned equation gives for $r\ne r^{\prime }$,

$\begin{array}{c}\hfill P^j\left[\frac{1}{h^2}\sum _{\eta \in {\mathcal{A}}_{|j}}{a}_{\eta |j}^{\left(2\right)}{G}_{\lambda |{\theta }^{\prime }}^{j+\eta ,k}\right]+\sum _{\mu \in \mathbb{Z}}{Q}_{\mu }^j\left[\frac{1}{h}\sum _{\eta \in {\mathcal{A}}_{|j}}{a}_{\eta |j}^{\left(1\right)}{G}_{\lambda -\mu |{\theta }^{\prime }}^{j+\eta ,k}\right]+\sum _{\mu \in \mathbb{Z}}{R}_{\lambda ,\mu }^j{G}_{\lambda -\mu |{\theta }^{\prime }}^{j,k}+O\left(h^{\xi }\right),\end{array}$

with $\xi$ the order of approximation, ${\mathcal{A}}_{|j}$ the set of neighbor site indices, and $a_{\eta |j}^{\left(n\right)}$ the respective coefficient (namely the finite difference coefficient included into a matrix representation ${\mathbf{A}}^{\left(n\right)}$—see Appendix B) of the $\eta$-th neighbor required to compute the n-th derivative up to a predetermined order of accuracy [42]; in the three-point stencil case, $\xi =2$. Both sets of neighbor indices and coefficients depend on the information of the site j under inspection; if, for example, we are at or near an interface, boundary, or discontinuity, then the strategy for choosing neighbors may differ; we might be interested in computing derivatives using only points in regions where it makes sense.

With some reorganization, the generated discrete DE can be regarded as a matrix multiplication. To see this, first we realize that by understanding $G_{\lambda |{\theta }^{\prime }}^{j,k}$ as the $(j,k)$-th element of a constructed matrix ${\mathbf{G}}_{\lambda |{\theta }^{\prime }}$—of size $N+1\times N+1$—we can envision a column matrix vector ${\mathbb{G}}_{{\theta }^{\prime }}$ that contains all $\lambda$-modes, or all of $\{{\mathbf{G}}_{\lambda |{\theta }^{\prime }}\forall \lambda \in \mathbb{Z}\}$, where all operations from the previous complex array equation are condensed into an equally conceived matrix $\mathbb{U}$ multiplying ${\mathbb{G}}_{{\theta }^{\prime }}$. The following is a view of ${\mathbb{G}}_{{\theta }^{\prime }}$,

(32)$\begin{array}{cc}\hfill & \hfill {\mathbb{G}}_{{\theta }^{\prime }}=\left(\begin{array}{c}⋮\\ {\mathbf{G}}_{-L|{\theta }^{\prime }}\\ ⋮\\ {\mathbf{G}}_{-1|{\theta }^{\prime }}\\ {\mathbf{G}}_{0|{\theta }^{\prime }}\\ {\mathbf{G}}_{1|{\theta }^{\prime }}\\ ⋮\\ {\mathbf{G}}_{L|{\theta }^{\prime }}\\ ⋮\end{array}\right),\mathrm{with}{\mathbf{G}}_{\lambda |{\theta }^{\prime }}=\left(\begin{array}{ccccc}{G}_{\lambda |{\theta }^{\prime }}^{0,0}& G_{\lambda |{\theta }^{\prime }}^{0,1}& \cdots & G_{\lambda |{\theta }^{\prime }}^{0,N-1}& G_{\lambda |{\theta }^{\prime }}^{0,N}\\ G_{\lambda |{\theta }^{\prime }}^{1,0}& G_{\lambda |{\theta }^{\prime }}^{1,1}& \cdots & G_{\lambda |{\theta }^{\prime }}^{1,N-1}& G_{\lambda |{\theta }^{\prime }}^{1,N}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ G_{\lambda |{\theta }^{\prime }}^{N-1,0}& G_{\lambda |{\theta }^{\prime }}^{N-1,1}& \cdots & G_{\lambda |{\theta }^{\prime }}^{N-1,N-1}& G_{\lambda |{\theta }^{\prime }}^{N-1,N}\\ G_{\lambda |{\theta }^{\prime }}^{N,0}& G_{\lambda |{\theta }^{\prime }}^{N,1}& \cdots & G_{\lambda |{\theta }^{\prime }}^{N,N-1}& G_{\lambda |{\theta }^{\prime }}^{N,N}\end{array}\right).\hfill \end{array}$

In principle, both matrices are infinitely large but for practical terms they will be truncated on both $\lambda$- and $\mu$-modes, as mentioned in the previous section. Despite this numerical simplification that will be carried out in the numerical analysis, the infinite matrix $\mathbb{U}$ has a well-defined structure, as will be detailed in Section 3.2.

Finally, the terms to the right of Equation (23) vanish for all $r\ne r^{\prime }$ leading us to believe that if $\mathbb{U}$ is invertible then the solution to the discrete Green function ${\mathbb{G}}_{{\theta }^{\prime }}$ is identically zero. However, attention should be paid at $r=r^{\prime }$ for its effect discards the trivial solution. Along with the other geometrical boundary conditions the problem will now have a unique solution. These boundary conditions will be addressed in Section 3.3.

3.2. Infinite Matrix $\mathbb{U}$

To understand the structure of $\mathbb{U}$ we turn to the set of operations for a particular $\lambda$-mode. Seeing as $\mathbb{U}$ is infinite we may encode rows by the integer value of the mode being solved and columns by the value of the mode being correlated. Thus taking row $\lambda$ from $\mathbb{U}$,

(33)$\begin{array}{cc}\hfill {\mathbb{U}}_{\lambda }{\mathbb{G}}_{{\theta }^{\prime }}& =\left(\underset{\cdots ,M,\cdots ,3,2,1⟵\mu |}{\cdots ,\stackrel{\mathrm{column}\lambda -\mu }{\overbrace{\frac{1}{h}{\mathbf{Q}}_{\mu }+{\mathbf{R}}_{\lambda ,\mu }}},\cdots ,}\underset{\mu =0}{\stackrel{\mathrm{column}\lambda }{\overbrace{\frac{1}{h^2}\mathbf{P}+\frac{1}{h}{\mathbf{Q}}_0+{\mathbf{R}}_{\lambda ,0}}}}\underset{|\mu ⟶-1,-2,\cdots ,-M,\cdots }{,\cdots ,\stackrel{\mathrm{column}\lambda -\mu }{\overbrace{\frac{1}{h}{\mathbf{Q}}_{\mu }+{\mathbf{R}}_{\lambda ,\mu }}},\cdots }\right)\hfill \\ \hfill & \phantom{=}\left(\begin{array}{cc}⋮& \\ {\mathbf{G}}_{\lambda -\mu |{\theta }^{\prime }}& \hfill \}\mathrm{row}\lambda -\mu \hfill \\ ⋮& \\ {\mathbf{G}}_{\lambda |{\theta }^{\prime }}& \hfill \}\mathrm{row}\lambda \hfill \\ ⋮& \\ {\mathbf{G}}_{\lambda -\mu |{\theta }^{\prime }}& \hfill \}\mathrm{row}\lambda -\mu \hfill \\ ⋮\end{array}\right),\hfill \end{array}$

with the following definitions for matrices $\mathbf{P}$, ${\mathbf{Q}}_{\mu }$, ${\mathbf{R}}_{\lambda ,\mu }$,

(34)$\begin{array}{cc}\hfill \mathbf{P}& \hfill =\left(\begin{array}{cccc}{P}^0& 0& \cdots & 0\\ 0& P^1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & P^N\end{array}\right)\times {\mathbf{A}}^{\left(2\right)},\hfill \end{array}$

(35)$\begin{array}{cc}\hfill {\mathbf{Q}}_{\mu }& \hfill =\left(\begin{array}{cccc}{Q}_{\mu }^0& 0& \cdots & 0\\ 0& Q_{\mu }^1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & Q_{\mu }^N\end{array}\right)\times {\mathbf{A}}^{\left(1\right)},\hfill \end{array}$

(36)$\begin{array}{cc}\hfill {\mathbf{R}}_{\lambda ,\mu }& \hfill =\left(\begin{array}{cccc}{R}_{\lambda ,\mu }^0& 0& \cdots & 0\\ 0& R_{\lambda ,\mu }^1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & R_{\lambda ,\mu }^N\end{array}\right).\hfill \end{array}$

One last remark on matrix $\mathbb{U}$ is that the density of non-zero entries is at most $3/N$ for the three-point stencil. For sufficiently large N, it will become essential to find a way to manage such sparsity for all speed-ups, data-compression, and efficiency in memory footprint.

3.3. Discrete Boundary Conditions

Using conditions detailed thoroughly in Section 2.2 and at end of Section 2.4 we are now in capacity of parameterizing the values of $G_{\lambda |{\theta }^{\prime }}^{j,k}$. This parametrization should further reflect the behavior of the $\delta$-function. The following are the conditions for the 3-point stencil: (i.) at $r=R_{\mathrm{ext}}$,

(37)$\begin{array}{cc}\hfill G_{\lambda |{\theta }^{\prime }}^{N,k}& \hfill =0\mathrm{for} \mathrm{DBC};\hfill \end{array}$

(38)$\begin{array}{cc}\hfill G_{\lambda |{\theta }^{\prime }}^{N+1,k}-G_{\lambda |{\theta }^{\prime }}^{N-1,k}& \hfill =0\mathrm{for} \mathrm{NBC},\hfill \end{array}$

(ii.) at $R_{\mathrm{int}}$ for the annulus,

(39)$\begin{array}{cc}\hfill G_{\lambda |{\theta }^{\prime }}^{0,k}& \hfill =0\mathrm{for} \mathrm{DBC};\hfill \end{array}$

(40)$\begin{array}{cc}\hfill G_{\lambda |{\theta }^{\prime }}^{1,k}-G_{\lambda |{\theta }^{\prime }}^{-1,k}& \hfill =0\mathrm{for} \mathrm{NBC},\hfill \end{array}$

(iii.) for the disc at $R_{\mathrm{int}}$ (disregarding $G_{\lambda }^{\prime }=0$ for now),

(41)$\begin{array}{cc}\hfill G_{0|{\theta }^{\prime }}^{1,k}-G_{0|{\theta }^{\prime }}^{-1,k}=0& \hfill \lambda =0,\hfill \end{array}$

(42)$\begin{array}{cc}\hfill G_{\lambda |{\theta }^{\prime }}^{0,k}=0& \hfill \lambda \ne 0,\hfill \end{array}$

and, finally, (iv.) at the interface $r=r^{\prime }$ the condition reads,

(43)$\begin{array}{c}\hfill \left(G_{\lambda |{\theta }^{\prime }}^{k_{>}+1,k}-G_{\lambda |{\theta }^{\prime }}^{k_{>}-1,k}\right)-\left(G_{\lambda |{\theta }^{\prime }}^{k_{<}+1,k}-G_{\lambda |{\theta }^{\prime }}^{k_{<}-1,k}\right)=\frac{2h}{r^k}{e}^{-i\lambda {\theta }^{\prime }}.\end{array}$

Here we have adopted the subscript convention of $<,>$ to refer to points to the left and right of the site of derivative evaluation. Note how all equations above reference and highlight a few fictitious points. The mesh points that lay outside or beyond the valid grid are $G_{\lambda |{\theta }^{\prime }}^{N+1,k}$, $G_{\lambda |{\theta }^{\prime }}^{-1,k}$, $G_{\lambda |{\theta }^{\prime }}^{k_{<}+1,k}$, and $G_{\lambda |{\theta }^{\prime }}^{k_{>}-1,k}$. These spurious terms must be dealt with and simplified in order to be able to incorporate readily all conditions.

3.4. A Non-Trivial Matrix Equation and a Solution

We will now show the explicit matrix equation associated 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 procedure for the three-point stencil and further discuss how to generalize for higher orders of approximation.

With the boundary relationships in mind, here in Equations (37), (39), (41) and (43), Equation (23) (multiplied by $-h^2$) equates partially to zero (when $r\ne r^{\prime }$) as,

$\begin{array}{cc}\hfill & \hfill 2P^j{G}_{\lambda |{\theta }^{\prime }}^{j,k}-P^j(G_{\lambda |{\theta }^{\prime }}^{j+1,k}+G_{\lambda |{\theta }^{\prime }}^{j-1,k})-h^2\sum _{\mu }{R}_{\lambda ,\mu }^j{G}_{\lambda -\mu |{\theta }^{\prime }}^{j,k}-\frac{h}{2}\sum _{\mu }{Q}_{\mu }^j(G_{\lambda -\mu |{\theta }^{\prime }}^{j+1,k}-G_{\lambda -\mu |{\theta }^{\prime }}^{j-1,k})=0,\hfill \end{array}$

where via Equation (43) the latter can be used to simplify both spurious terms (appearing at $r=r^{\prime }$) $G_{\lambda |{\theta }^{\prime }}^{k_{<}+1,k}$ and $G_{\lambda |{\theta }^{\prime }}^{k_{>}-1,k}$. After crossing out these terms by iterative substitution we obtain a generalized expression for the above valid for almost every point in the grid. The general discrete equation yields for $r^{\prime }>0$,

(44)$\begin{array}{cc}\hfill & \hfill 2P^j{G}_{\lambda |{\theta }^{\prime }}^{j,k}-P^j(G_{\lambda |{\theta }^{\prime }}^{j+1,k}+G_{\lambda |{\theta }^{\prime }}^{j-1,k})-\frac{h}{2}\sum _{\mu }{Q}_{\mu }^j(G_{\lambda -\mu |{\theta }^{\prime }}^{j+1,k}-G_{\lambda -\mu |{\theta }^{\prime }}^{j-1,k})-h^2\sum _{\mu }{R}_{\lambda ,\mu }^j{G}_{\lambda -\mu |{\theta }^{\prime }}^{j,k}\hfill \\ \hfill & \hfill =-h{r}^j{\delta }^{j,k}{e}^{-i\lambda {\theta }^{\prime }}\left\{1-\frac{h^2}{4{\left(r^k\right)}^2}{[1+r^k{f}_r(r^k,{\theta }^{\prime })]}^2\left[1-{\delta }_{\lambda ,0}\right]\right\},\hfill \end{array}$

where the new term that accounts for the boundary condition at $r=r^{\prime }$ has appeared. Due to the absence of a left-hand limit as $r^{\prime }=0$, according to Equation (43), this term is exactly $-h{r}^j{\delta }^{j,k}{e}^{-i\lambda {\theta }^{\prime }}$ at the origin. Actually, this condition holds for the mode $\lambda =0$ in general due to translational invariance—this invariance is clearly absent for the other modes. The terms composing the right-hand side of last equation can be viewed as of order of mesh-size or order of radial distance from the origin as follows,

• 1.. $-\frac{h^3}{4}{f}_r(r^k,{\theta }^{\prime })$, a surprising third order correction due to the vector field appearing after substituting the interface difference in derivatives.

• 2.. $h{r}^k=h{R}_{\mathrm{int}}+h^{2}k$, the leading order that substituted yields a first order constant term and a second order increasing term.

• 3.. $-\frac{h^3}{4}\frac{1}{r^k}=-\frac{h^3}{4}\frac{1}{R_{\mathrm{int}}+hk}$, a negative term significant closer to the origin. As expected, the behavior of the discrete version near zero validates our previous choice of boundary condition for the disc.

• 4.. For $r^{\prime }=0$, we must implement a cut-off such that $r^0=ϵ>0$ instead of zero to avoid numerical divergences. The choice for $ϵ$ will be discussed below.

Because the error in the differential equation is of $O\left(h^4\right)$, we should incorporate all terms to the calculation. However, we will neglect the higher order term—first term—since this will simplify our calculations of $\psi \left(\mathbf{r}\right)$.

This final expression is valid everywhere including the controversial $j=0,N$ points, where either Dirichlet or Neumann conditions complete Equation (44) at the borders. In those two cases, substitutions must take place following Equations (37), (39) and (41). After replacements, and due to the nature of derivative calculation in the three-point stencil, rows from $\mathbb{U}$ corresponding to exterior and interior borders are modified. See the substitution rules in Table 2 and Table 3.

We now define our complete matrix system as $\mathbb{U}·{\mathbb{G}}_{{\theta }^{\prime }}=-h\mathbb{V}·{\mathbb{E}}_{{\theta }^{\prime }}$. Among other things, the right-hand side accounts for the contribution of Dirac’s distribution. The two additional definitions appearing correspond to first a distance parameter generalized into ${\alpha }_{\lambda }^j$, a new object that incorporates the boundary conditions at both $j=0$ and $j=N$. Notice, for instance, that keeping the term $r^j$ at every point does not explain the vanishing of the Green function at the boundaries when DBC are considered, neither does it describe the correct behavior at $r=0$ for a disk. Actually, when the last condition is considered, an ultraviolet cut-off $ϵ$—such that $ϵ\to 0$—must be introduced to avoid divergences, as seen in previous works [14,15,16]. Such a cut-off is not surprising, as the 2D Green distribution has a natural divergence at $\mathbf{r}={\mathbf{r}}^{\prime }$ and a logarithmic behavior near the origin when ${\mathbf{r}}^{\prime }=0$. Although the appearance of this divergence can easily be visualized after studying the behavior of Equation (44) at $j=0$ for the mode $\lambda =0$ in a disk, its existence at any point—also for an annulus—is guaranteed by the infinite number of $\lambda$-modes that must be summed up to obtain an exact solution. Therefore, it is not surprising that $ϵ$ and h are related—see Section 4 for more details.

Matrix terms are written as,

(45)$\begin{array}{cc}\hfill {\mathbb{V}}_{\lambda ,\mu }^{j,k}& \hfill ={\alpha }_{\lambda }^j{\delta }_{\lambda ,\mu }{\delta }^{j,k},\hfill \end{array}$

(46)$\begin{array}{cc}\hfill {\mathbb{E}}_{\lambda |{\theta }^{\prime }}^{j,k}& \hfill =e^{-i\lambda {\theta }^{\prime }}{\delta }^{j,k},\hfill \end{array}$

and the solution to the $\lambda$-modes matrix ${\mathbb{G}}_{{\theta }^{\prime }}$ is,

(47)$\begin{array}{c}\hfill {\mathbb{G}}_{{\theta }^{\prime }}=-h\mathbb{A}·{\mathbb{E}}_{{\theta }^{\prime }},\end{array}$

where we have defined $\mathbb{A}={\mathbb{U}}^{-1}·\mathbb{V}$ assuming that $\mathbb{U}$ is invertible.

Matrix $\mathbb{A}$ was declared because it has interesting symmetry properties that will be discussed in the next section. Table 2, Table 3, Table 4 and Table 5 outline how to fill the matrix elements of the objects we have described.

3.5. The Parameter $\alpha$ and the Symmetry of $\mathbb{A}$

A closed relation can be found for the matrix describing the entire Green function. Using the results from previous section, it is

(48)$\begin{array}{c}\hfill G_{\theta ,{\theta }^{\prime }}^{j,k}=-\frac{h}{2\pi }\sum _{\lambda ,\mu }{e}^{i\lambda \theta }{e}^{-i\mu {\theta }^{\prime }}{\mathbb{A}}_{\lambda ,\mu }^{j,k},\end{array}$

where ${\mathbb{A}}_{\lambda ,\mu }^{j,k}={\left[{\mathbb{U}}^{-1}\right]}_{\lambda ,\mu }^{j,k}{\alpha }_{\mu }^k$. As previously mentioned, the parameter ${\alpha }_{\lambda }^k$ generalizes the radial parameter $r^k$, including the boundary conditions. On the other hand, it is worthwhile to state the symmetry conditions that $\mathbb{A}$ satisfies (a more detailed derivation can be found in Appendix D; $\overline{z}$ denotes the complex conjugate of z.)

(49)$\begin{array}{cc}\hfill & \hfill Im\left({\mathbb{A}}_{0,0}^{j,k}\right)=0,\hfill \end{array}$

(50)$\begin{array}{cc}& \hfill {\mathbb{A}}_{-\lambda ,0}^{j,k}=\overline{{\mathbb{A}}_{\lambda ,0}^{j,k}},{\mathbb{A}}_{0,-\mu }^{j,k}=\overline{{\mathbb{A}}_{0,\mu }^{j,k}},\hfill \end{array}$

(51)$\begin{array}{cc}\hfill & \hfill {\mathbb{A}}_{-\lambda ,\mu }^{j,k}=\overline{{\mathbb{A}}_{\lambda ,-\mu }^{j,k}},{\mathbb{A}}_{-\lambda ,-\mu }^{j,k}=\overline{{\mathbb{A}}_{\lambda ,\mu }^{j,k}}.\hfill \end{array}$

Using Table 5, it is easy to see that the matrix elements ${\left[{\mathbb{U}}^{-1}\right]}_{\lambda ,\mu }^{j,k}$ satisfy the same symmetry properties.

3.6. The Algorithm

The algorithm for a numerical solution can be summarized as follows:

• 1.. The values of L and M are determined according to the required level of approximation.

• 2.. We fill all elements described in Table 1; the matrix elements ${\mathbb{U}}_{\lambda ,\mu }^{j,k}$ are filled by blocks using the rules shown in Table 2, Table 3 and Table 4. Matrix elements ${\alpha }_{\lambda }^j$ are also filled according to Table 5.

• 3.. Matrix $\mathbb{U}$ is inverted and so matrix $\mathbb{A}$ is computed.

• 4.. The Green function is computed according to Equation (48).

Using previous results, we can deduce a closed form for $\psi \left(\mathbf{r}\right)$ for both DBC and NBC using the conventions stated in Equations (8) and (14). Although there are many ways to perform the integrals stated in previous equations, and the reader can choose the method that he or she prefers, a sketch of these solutions, using the trapezoid rule, is shown in Appendix F and Appendix G.

3.7. Particular Cases and Properties

We will analyze some particular cases that can be deduced from the procedure explained above. We will start focusing on the one-dimensional case.

• One-Dimensional Case

The analysis of a Green function in one dimension requires an appropriate definition of a general DE obeyed by the Green function $G(x,x^{\prime })$. Unfortunately, a direct analysis of the results by studying Equation (23) is not straightforward due to the clear differences between the Laplacians in Cartesian and polar coordinates. Let us imagine a general second order DE of the form ${\mathcal{L}}_x\psi \left(x\right)=\varphi \left(x\right)$, where the Green function satisfies the relation

(52)$\begin{array}{cc}\hfill {\mathcal{L}}_{x}G(x,x^{\prime })& \hfill =P\left(x\right)\frac{d^{2}G}{d{x}^2}+Q\left(x\right)\frac{dG}{dx}+R\left(x\right)G=\beta \delta (x-x^{\prime }).\hfill \end{array}$

The function $P\left(x\right)$ might not be necessary, as it can be eliminated by division, but its inclusion allows us to have a more general analysis. The constant therm $\beta$ seems clumsily placed, as its value is usually 1. Nonetheless, some formalisms define the Green function by means of the operator ${\mathcal{L}}_{x}G(x,x^{\prime })=-\delta (x-x^{\prime })$, thus introducing a 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 for Equation (52), which is

(53)$\begin{array}{cc}\hfill & \hfill (2P^j-h^2{R}^j)G^{j,k}-(P^j+\frac{h}{2}{Q}^j)G^{j+1,k}-(P^j-\frac{h}{2}{Q}^j)G^{j-1,k}=-h\beta {\delta }^{j,k}.\hfill \end{array}$

Having confined the system within the domain $x\in [x^0,x^N]$, our step size is now $h=(x^N-x^0)/N$.

From this point on, we can apply the results obtained for the two-dimensional problem in this study. Notice that, in the absence of modes that account for the angular dependence, we can always say that ${\mathbb{A}}_{\lambda ,\mu }^{j,k}={\mathbb{A}}_{\lambda ,\mu }^{j,k}{\delta }_{\lambda ,0}{\delta }_{\mu ,0}$. Therefore, Equation (48) becomes

(54)$\begin{array}{c}\hfill G^{jk}=-h{\mathbb{A}}^{j,k},\mathrm{where}{\mathbb{A}}^{j,k}={\left[{\mathbb{U}}^{-1}\right]}^{j,k}{\alpha }^k\end{array}$

and ${\alpha }^k$ accounts for the boundary conditions. By making the association $x^j=x^0+hj$, the elements described in Equation (54), for $1\le j\le N-1$, are now filled using the following rules:

• 1.. ${\mathbb{U}}^{j,j}=2P^j-h^2{R}^j$.

• 2.. ${\mathbb{U}}^{j,j\pm 1}=-P^j\mp \frac{h}{2}{Q}^j$.

• 3.. ${\mathbb{U}}^{0,0}={\mathbb{U}}^{N,N}=1$ for DBC. For NBC: ${\mathbb{U}}^{0,0}=2q^0-h^2{b}^0$, ${\mathbb{U}}^{N,N}=2P^N-h^2{R}^N$, ${\mathbb{U}}^{0,1}=-2P^1$ and ${\mathbb{U}}^{N,N-1}=-2P^{N-1}$.

• 4.. ${\alpha }^j=\beta$.

• 5.. ${\alpha }^0={\alpha }^N=0$ for DBC. For NBC: ${\alpha }^0=\beta$ and ${\alpha }^N=\beta$.

The weight function, which guarantees the symmetry condition $G^{kj}=w^{jk}{G}^{jk}$ takes the form

(55)$\begin{array}{c}\hfill w(x^{\prime },x)=\frac{e^{U\left(x^{\prime }\right)}}{e^{U\left(x\right)}},U\left(z\right)=\frac{1}{P\left(z\right)}exp\left[{\int }_{z_0}^z\frac{Q\left(y\right)}{P\left(y\right)}dy\right],\end{array}$

where $z_0$ is an irrelevant constant. Solutions for $\psi \left(x\right)$ with both DBC and NBC using the trapezoid rule as method of integration are shown in Appendix G.

• Monopole-like Case

This takes place when both $\overrightarrow{f}\left(\mathbf{r}\right)$ and $g\left(\mathbf{r}\right)$ have no significant angular dependence, so the mode $\mu =0$ is their only relevant contribution; this implies that $Q_{\mu }^j=R_{\lambda ,\mu }^j=0$ for $\mu \ne 0$. The off-diagonal matrices ${\mathbb{U}}_{\lambda ,\lambda -\mu }^{j,k}$ thus vanish—this leads to a block diagonal $\mathbb{U}$ matrix—and so the system becomes separable in the radial and angular variables. Having now the relation ${\mathbb{A}}_{\lambda ,\mu }^{j,k}={\mathbb{A}}_{\lambda ,\mu }^{j,k}{\delta }_{\lambda ,\mu }$, Equation (48) reduces to

(56)$\begin{array}{c}\hfill G_{\theta {\theta }^{\prime }}^{jk}=-\frac{h}{2\pi }\sum _{\lambda }{e}^{i\lambda (\theta -{\theta }^{\prime })}{\mathbb{A}}_{\lambda ,\lambda }^{j,k}.\end{array}$

Notice that each mode can now be solved independently.

3.8. Beyond the Three-Point Stencil

As mentioned, we only showed an explicit analysis for a three-point stencil approximation. This method can be generalized to include the contribution of more neighbors in the derivative terms, i.e., higher order stencils that provide more accurate degrees of approximation in h. In spite of its simplicity, the three-point stencil has the great advantage that the spurious terms that arise from the boundary conditions can be eliminated in a simple fashion.

The description of the system with a five-point stencil, for instance, will increase the amount of terms different from zero in $\mathbb{U}$—for example, terms of form ${\mathbb{U}}_{\lambda ,\lambda -\mu }^{j,j\pm 2}$ will provide non-zero contributions. Having a higher degree of approximation, that demands the inclusion of more non-trivial matrix terms, the grid size can be reduced. Although there is no guarantee that the inversion process is optimized in time when the contributions of more neighbors are included, as the matrices are highly sparse, there is a clear optimization of memory storage.

Yet a great disadvantage that higher stencils inherit is the elimination of the spurious terms that come from the boundary conditions. For instance, when we deal with the condition at $j=k$ ($r=r^{\prime }$), Equation (43) will include more coefficients outside the grid, so the recurrence relation that is obtained will not be able to eliminate all of them—at least, using the same procedure we implemented. Therefore, a different approach must be performed. A possible solution could be expanding the derivatives around a point different from the center, so avoiding the spurious terms; this process is studied in detail in [42]. Nonetheless, this could be discussed in a future study.

4. Numerical Results

We will use the formalism described above to solve some particular examples. This section is particularly focused on presenting our results; such results will be discussed in the following section. Three examples will be presented.

• Example 1: a one-dimensional case

As a first example, let us study a one-dimensional system with a known analytical solution, useful to test the formalism we have described. Let us suppose we want to solve the DE in the domain $[0,4]$

(57)$\begin{array}{c}\hfill x^2{\psi }^{″}\left(x\right)+x{\psi }^{\prime }\left(x\right)+(x^2-4)\psi \left(x\right)=J_4\left(x\right),\end{array}$

where $J_n\left(x\right)$ and $Y_n\left(x\right)$ are the Bessel functions of first and second kind of order n. Using Equation (55), we can easily deduce that $w(x^{\prime },x)=x/x^{\prime }$.

The conditions ${\psi }^0=0$ and ${\psi }^N=2$ (Dirichlet), lead to the analytical solution,

(58)$\begin{array}{cc}\hfill \psi {\left(x\right)}^{\mathrm{DBC}}=\frac{1}{24J_2\left(4\right)}[& \hfill J_2\left(x\right)\left(48-2J_4\left(4\right)++\pi x{J}_2\left(4\right)J_4\left(x\right)Y_1\left(x\right)\right)\hfill \\ \hfill & \hfill -\pi x{J}_2\left(4\right)J_1\left(x\right)J_4\left(x\right)Y_2\left(x\right)].\hfill \end{array}$

Conversely, with conditions ${f^{\prime }}^0=0$ and ${f^{\prime }}^N=2$ (Neumann), the analytical solution yields,

(59)$\begin{array}{cc}\hfill \psi {\left(x\right)}^{\mathrm{NBC}}=\frac{1}{C\left(x\right)}[& \hfill \left(2J_0\left(4\right)+3J_1\left(4\right)\right)\left(x(x^2-24)J_0\left(x\right)-8(x^2-6)J_1\left(x\right)\right)\hfill \\ \hfill & \hfill +x^3{J}_2\left(x\right)\left(96+2J_0\left(4\right)-3J_1\left(4\right)\right)],\hfill \end{array}$

with $C\left(x\right)=24x^3\left(J_1\left(4\right)-J_3\left(4\right)\right)$.

Analytical and numerical results are compared for both cases in Figure 2 and Table 6—see Equations (A25a) and (A25b) for explicit expressions using a numerical approach.

• Example 2: the two-dimensional Helmholtz equation with imaginary wave number

We now shift our attention to solve a two-dimensional system. Let us consider the DE

(60)$\begin{array}{c}\hfill ({\overrightarrow{\nabla }}^2-m^2)G(\mathbf{r},{\mathbf{r}}^{\prime })=\delta (\mathbf{r}-{\mathbf{r}}^{\prime }).\end{array}$

In the absence of the vector field $\overrightarrow{f}\left(\mathbf{r}\right)$, we conclude that $w(\mathbf{r},{\mathbf{r}}^{\prime })=1$; besides, the system is separable in the radial and angular coordinates. The solution to last equation confined in a large disc of radius $r^N=R_{\mathrm{ext}}$ with Dirichlet boundary conditions can be found analytically. For $\theta ={\theta }^{\prime }$ and different values of $r^{\prime }$, solutions are presented in Figure 3.

Last numerical result will now be used to solve a non-homogeneous equation of the form $({\overrightarrow{\nabla }}^2-m^2)\psi \left(\mathbf{r}\right)=\varphi \left(\mathbf{r}\right)$, whose general solution is provided in Appendix G with $r^0=R_{\mathrm{int}}=0$.

Now, let us consider $\varphi \left(\mathbf{r}\right)$ to be a function defined over a disc of radius $R=10$ and study the two following cases:

• (a). ${\psi }^{\left(a\right)}(r,\theta )$, with $\varphi \left(\mathbf{r}\right)=-\frac{1}{10}rsin\theta$, and $\psi (R,\theta )=2$.

• (b). ${\psi }^{\left(b\right)}(r,\theta )$, with $\varphi \left(\mathbf{r}\right)=\left\{\begin{array}{ccc}{r}^{-1/2}& ,& 0\le \theta <\pi \\ -r^{-1/2}& ,& \pi \le \theta <2\pi \end{array}\right.$, and

  $\psi (R,\theta )=\left\{\begin{array}{ccc}1& ,& 0\le \theta <\pi \\ -1& ,& \pi \le \theta <2\pi \end{array}\right.$.

Solutions to ${\psi }^{\left(a\right)}(r,\theta )$ and ${\psi }^{\left(b\right)}(r,\theta )$ for some angles are shown in Figure 4 and Figure 5, respectively.

It is important to mention that Equation (60) is a particular case of a DE that describes the behavior of a two-dimensional two-component plasma in thermal equilibrium at a fixed temperature in the absence of an external electric field. Here the function $\psi \left(\mathbf{r}\right)$ is naturally identified with $G(\mathbf{r},{\mathbf{r}}^{\prime })$ (actually, a Green function), which generates the n-body correlation functions of the system. When the plasma is confined in a disc (or annulus) and it is subjected to an external field, $\overrightarrow{E}\left(\mathbf{r}\right)$, the system can be solved numerically with the presented method, once the identifications $g\left(\mathbf{r}\right)=-m^2-E^2\left(\mathbf{r}\right)$, $f_x\left(\mathbf{r}\right)=-2i{E}_y\left(\mathbf{r}\right)$, and $f_y\left(\mathbf{r}\right)=2i{E}_x\left(\mathbf{r}\right)$ are performed; more details can be found in [14,15,16].

• Example 3: a pedagogical example

Now let us apply the same formalism to solve another two-dimensional problem. Let us suppose that we want to find the Green function associated with the two-dimensional DE $[\overrightarrow{\nabla }+\overrightarrow{\nabla }Z\left(\mathbf{r}\right)]·\overrightarrow{\nabla }\psi \left(\mathbf{r}\right)=\varphi \left(\mathbf{r}\right)$, where $Z(x,y)=2x^2{y}^2$. After transforming the system to polar coordinates, we can see that $f_r=r^3(1-cos4\theta )$ and $f_{\theta }=r^{3}sin4\theta$. The elements defined in Table 1 now become

(61a)$\begin{array}{cc}\hfill Q_{\mu }^j& \hfill =r^j{\delta }_{0,\mu }+{\left(r^j\right)}^5\left[{\delta }_{0,\mu }-\frac{1}{2}({\delta }_{4,\mu }+{\delta }_{-4,\mu })\right],\hfill \end{array}$

(61b)$\begin{array}{cc}\hfill R_{\lambda ,\mu }^j& \hfill =-{\lambda }^2{\delta }_{0,\mu }+\frac{1}{2}{\left(r^j\right)}^4(\lambda -\mu )({\delta }_{4,\mu }-{\delta }_{-4,\mu }).\hfill \end{array}$

The system will be confined in an annulus or internal radius $R_{\mathrm{int}}=1$ and external radius $R_{\mathrm{ext}}=2$. Figure 6 shows the Green function for DBC and some particular parameters but different values of L. Figure 7 shows the results for different parameters under NBC.

5. Discussion

It is worthwhile to discuss the results found in the previous section. We will analyze the three examples presented above in respective order followed by an explanation on how the method can be adapted to solve the diffusion equation.

Example 1.

In the one-dimensional case it is interesting to contrast the numerical solutions shown in Figure 2 using the weight function and the one that arises without the weight function formalism—performing the replacement $w^{jk}{G}^{jk}\to G^{kj}$ (graphs that illustrate the results using the weight function were not shown, but comparisons will be given. Interestingly, from Table 6, we conclude that the introduction of the weight function leads to a more accurate result in this case. Notice that the use of a modest grid sizes ($32\le N\le 512$) leads to very accurate results.

Example 2.

This two-dimensional example has an analytical solution. Adapting the result found in [16], we deduce that the Green function associated with Equation (60) is

(62) $\begin{array}{c}\hfill G(\mathbf{r},{\mathbf{r}}^{\prime })=-\frac{m^2}{2\pi }\sum _{\lambda }{e}^{i\lambda (\theta -{\theta }^{\prime })}\left[I_{\lambda }\left(m{r}_{<}\right)K_{\lambda }\left(m{r}_{>}\right)-t_{\lambda }\left(R\right)I_{\lambda }\left(mr\right)I_{\lambda }\left(m{r}^{\prime }\right)\right],\end{array}$

where $r_{>}$ and $r_{<}$ are the maximum and minimum between r and $r^{\prime }$, $I_{\lambda }\left(x\right)$ and $K_{\lambda }\left(x\right)$ are the well-known modified Bessel function of the second kind and $t_{\lambda }\left(R\right)=\frac{K_{\lambda }\left(mR\right)}{I_{\lambda }\left(mR\right)}$. Taking a look at Equation (62), we deduce that the Green function diverges—many distributions are formally infinite. Actually, the first term of previous sum can be reduced to $-\frac{m^2}{2\pi }{K}_0\left(m\right|\mathbf{r}-{\mathbf{r}}^{\prime }\left|\right)$.

Notice that the Green function diverges logarithmically, in the limit $\mathbf{r}\to {\mathbf{r}}^{\prime }$, as $K_0{\left(x\right)}_{x\to 0}\sim ln(2/x)-\gamma$, with $\gamma$ the Euler Mascheroni constant. A cut-off s, which represents a minimum separation distance between $\mathbf{r}$ and ${\mathbf{r}}^{\prime }$ [14], is usually introduced to address this divergence. In turn, s might be related to L and $ϵ$.

We now implement the numerical analysis to verify last solution noticing that $P^j={\left(r^j\right)}^2$, $Q_{\mu }^j\left(r\right)={\alpha }^j=r^j{\delta }_{\mu ,0}$, and $R_{\lambda ,\mu }^j=-({\lambda }^2+m^2{P}^j){\delta }_{\mu ,0}$.

The following step is determining an appropriate value for $ϵ$. From Figure 3, we see that the solution for a fixed value of $r^{\prime }$ shows a peak at $\mathbf{r}={\mathbf{r}}^{\prime }$. As we sum up all modes the magnitude of the height of the peaks must be infinite. However, the introduction of L guarantees the peaks to be finite. The height of the peak at $r=r^{\prime }=0$ is associated with $ϵ$, as the functions $K_{\lambda }\left(r\right)$ diverge at $r=0$. The cut-off $ϵ$ is chosen in such way that the height of the peaks in the neighborhood of $r=r^{\prime }=0$—i.e., $|r-r^{\prime }|=O\left(h\right)$—are close enough. We found that for $N\ge 2^7$, $ϵ\simeq 0.25h$. Surprisingly, we found that $ϵ$ does not depend on L for large enough L.

We recall that this is an approximation; the exact solution for the distribution is found in the limits $ϵ\to 0$ and $L\to \infty$. In addition to the graphical comparison between the numerical and analytical results already shown in Figure 3, Table 7 also presents numerical details that ratify the accuracy of our results.

Regarding the inhomogeneous solutions shown in Figure 4 and Figure 5, we do not have an analytical solution to test the algorithm; however, results show to be consistent. Particularly, the function ${\psi }^{\left(a\right)}(r,\theta )$ shows maxima for similar values of r—the value of such maxima gradually decreases with the angular direction; interestingly, the location of the critical points seems to occur for approximately the same value of r. Since $\varphi \left(\mathbf{r}\right)$ obeys a diffusion equation, this suggests that the concentration of the analyzed quantity shows maximum (or minimum) values at a certain radial location; this critical concentration, nonetheless, is favored for larger angles—this is consistent because of the non-homogeneous function $\varphi \left(\mathbf{r}\right)=-\frac{1}{10}rsin\theta$ that has been chosen. As expected, ${\psi }^{\left(a\right)}(r,\theta )$ take the required values at the borders for every angle. On the other hand, ${\psi }^{\left(b\right)}(r,\theta )$ presents the same solution for two sets of angles: $0\le \theta <\pi$ and $\pi \le \theta <2\pi$, which is expected due to the form chosen for $\varphi \left(\mathbf{r}\right)$—notice that the only difference between such two sets of conditions is a global sign; this difference is clearly reflected in Figure 5. Interestingly, $\varphi \left(\mathbf{r}\right)$ is finite at $r=0$, in spite of the clear divergence of $\varphi \left(\mathbf{r}\right)$ at this value. However, such difference is weak enough as it goes like $r^{-1/2}$. Like in previous example, the presence of a critical point is evident.

Example 3.

Notice from Figure 6 how the Green functions become zero at the borders and their derivatives present a discontinuity at $r=r^{\prime }$. As expected, the larger L, the larger the magnitude of the value at that point; however, this value decreases as $|\theta -{\theta }^{\prime }|$ increases; under this condition, the behavior of the Green function is smoother, but it seems to be too much more sensitive respect to the value of L. Interestingly, for the parameters that have been analyzed, the Green function only shows a “mirror symmetry” respect to $r^{\prime }$ as $\theta ={\theta }^{\prime }=0$. For NBC, as illustrated in Figure 7, something similar happens. However, the derivatives at the borders are now the ones which tend no be zero. Although this is clear as $r\to 2$, the asymptotic behavior toward zero close to the inner border can be appreciated. The values of the Green function at the border, however, clearly differ for different parameters. Similarly to the case where DBC where analyzed, the “mirror symmetry” around $r^{\prime }$ is also absent (also when $\theta ={\theta }^{\prime }=0$); the variation of the Green function is also larger as $\theta \simeq {\theta }^{\prime }$, showing sharper peaks.

When using the Green function to solve a particular inhomogeneous equation, it is clear that $\overrightarrow{\nabla }\times f=0$, so the weight function exists. It is

(63)$\begin{array}{c}\hfill w({\mathbf{r}}^{\prime },\mathbf{r})=\frac{e^{\frac{1}{4}{r^{\prime }}^4(1-cos4{\theta }^{\prime })}}{e^{\frac{1}{4}{r}^4(1-cos4\theta )}}=\frac{{\sum }_{\lambda }{w}_{\lambda }\left(r^{\prime }\right)e^{i\lambda {\theta }^{\prime }}}{w\left(\mathbf{r}\right)},\end{array}$

whose only non-vanishing modes in discrete coordinates are given by the relation $w_{4\lambda }^j={(-1)}^{\lambda }{e}^{\frac{1}{4}{\left(r^j\right)}^4}{I}_{\lambda }\left(\frac{1}{4}{\left(r^j\right)}^4\right)$. Although we are not interested in finding $\psi \left(\mathbf{r}\right)$ for a particular boundary problem, all the steps are carried out to accomplish this goal.

• The Stationary Diffusion Equation

It is worthwhile to describe how our formalism can be adapted to solve the diffusion equation at “thermal” equilibrium. Let $\psi \left(\mathbf{r}\right)$ and $D\left(\mathbf{r}\right)$ represent the density of the diffusion material and the anisotropic diffusion coefficient, respectively. In the stationary regime, $\psi \left(\mathbf{r}\right)$ satisfies the DE

(64)$\begin{array}{c}\hfill D\left(\mathbf{r}\right){\overrightarrow{\nabla }}^2\psi \left(\mathbf{r}\right)+\overrightarrow{\nabla }D\left(\mathbf{r}\right)·\overrightarrow{\nabla }\psi \left(\mathbf{r}\right)=0.\end{array}$

Although Equation (64) does not have the standard form shown in Equation (2), after dividing Equation (64) by $D\left(\mathbf{r}\right)$ and defining $\overrightarrow{f}\left(\mathbf{r}\right)$ as $\overrightarrow{f}\left(\mathbf{r}\right)=\frac{1}{D\left(\mathbf{r}\right)}\overrightarrow{\nabla }D\left(\mathbf{r}\right)$, the standard form can be obtained. (The diffusion coefficient is assumed to be well-behaved within the annulus or disc. However, $\overrightarrow{f}$ can have poles within the same domain.) The weight function is now guaranteed to exist, as $\overrightarrow{\nabla }\times \overrightarrow{f}=-\frac{1}{D^2}\overrightarrow{\nabla }D\times \overrightarrow{\nabla }D+\frac{1}{D}\overrightarrow{\nabla }\times \overrightarrow{\nabla }D=0$. Actually, $w(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{D\left(\mathbf{r}\right)}{D\left({\mathbf{r}}^{\prime }\right)}$.

The viability of our method to solve Equation (64) depends on the particular form of the diffusion coefficient. The following possibilities may arise: (a) $\overrightarrow{f}$ has no poles in the two dimensional domain; (b) $\overrightarrow{f}$ has a divergence in $r=0$ that can be eliminated once $\overrightarrow{f}$ is multiplied by r; (c) the divergence at $r=0$—or any other radial divergence—previously discussed still persists after multiplication by r; and (d) $D^{-1}$ has poles for some $\theta \in [0,2\pi )$.

The cases (a) and (b) can be solved with the regular procedure we have described; the modes $f_{r\mu }$ and $f_{\theta \mu }$ are well-behaved and so the elements $Q_{\mu }^j$ and $R_{\lambda ,\mu }^j$ defined in Table 1 exist. The possibility stated in (c) demands a redefinition of $\overrightarrow{f}$ to eliminate any possible radial divergence; however, this redefinition does not guarantee the existence of the weight function. The situation described in (d) is problematic, as some of the modes $f_{r\mu }$ and $f_{\theta \mu }$ are divergent. The last situation is alleviated by working with the original DE, Equation (64); nonetheless, the process that we must follow to solve a system whose mathematical form differs from Equation (2) has not been described in this work.

Similar analysis can be performed as we deal with the Poisson’s equation associated with electrostatic potential in an anisotropic media, among others.

6. Conclusions

In this paper we have analyzed the Green function formalism and studied under which conditions such mechanisms can be used to obtain the solution of an inhomogeneous DE. We consider that the main contribution of this research has been the proposition of a numeric algorithm to solve a family of partial differential equations for two-dimensional systems in polar coordinates, the Green function formalism was used for this purpose.

After decomposing the Green function as a sum of Fourier modes, an infinite set of coupled second order DE for the radial variable is found. Although such set decouples when the initial DE is separable in the radial and angular variables, the coupling in the modes arises as the vector field $\overrightarrow{f}\left(\mathbf{r}\right)$ and the scalar function $g\left(\mathbf{r}\right)$ are expressed as a sum of Fourier modes. The presented analysis allows us to state the following highlights: (a) we found that there exists a function, which we called the weight function, that makes the Liouville operator self-adjoint. This function also defines the symmetry properties of the Green function (how it is transformed under the exchange of $\mathbf{r}$ and ${\mathbf{r}}^{\prime }$) and is guaranteed to exist as $\overrightarrow{\nabla }\times \overrightarrow{f}=0$; (b) a non-separable DE in r and $\theta$ leads to a large matrix system that is non-diagonal by blocks; (c) there exist a natural diverge due to the usual properties of the Dirac Delta distribution; however, this divergence can by absorbed only in the $\lambda =0$-mode; and (d) the $\lambda$-modes satisfy DBC at $r=0$ for a disk geometry.

An algorithm to solve the Green function associated with a general class of Liouville operator was solved using a FEM. We used a simple three-point stencil approach to approximate the solution and focused on both Dirichlet and Neumann boundary conditions. A set of approximations was made, which included a truncation of the infinite number of modes, a minimum distance when the system is confined in a disc, and the discard of the term $-\frac{h^3}{4}{f}_r(r^k,{\theta }^{\prime })$. Although the first two approximations are well-justified because the Green function has a natural divergence, the last one was performed by convenience (anyway, it provides a very small contribution).

The algorithm was verified by comparing with known results and obtaining very small percentage errors. An additional example whose solution cannot be found by means of the regular algorithms was shown.

We consider that the presented method is a useful attempt to solve Green functions of operators whose radial and angular variables cannot be separated. However, we expect this algorithm can be improved by other authors in the future to obtain more accuracy without the need of creating huge matrix systems, which demand large storage memory and computational time. Some of the improvements may include the implementation of the method for higher order stencils, an optimized calculation either mathematically or numerically of $ϵ$, simplified formulas for the calculation of $\psi \left(\mathbf{r}\right)$ or “on the go” algorithms that do not require the inversion of the matrix or the storage of temporal information.

Author Contributions

A.F. was in charge of writing the article, developing the introduction and the main steps of the methodology, working on the numerical codes, and obtaining the results. Additionally, he proposed the project that did lead to the realization of this work. J.P.M. made relevant corrections and suggestions to the main document, proved some of the equations that are stated, help to consolidate the numerical codes, and gave some ideas to extend this work to future research projects. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript: Finite Elements Metyhod (FEM), Finite Difference Method (FDM), Diichlet Boundary Conditions (DBC), Neumann Boundary Conditions (NBC), Differential Equation (DE), Boundary Conditions (b.c.), Partial Differential Equations (PDE).

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Contours used in the path integration: to the left [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] to the right. The radius and thickness are chosen purposely as [Forumla omitted. See PDF.] to take the [Forumla omitted. See PDF.] limit.

Figure 2. Left: solution to the DE given by Equation (57) with the initial conditions [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. Right: solution to the DE given by Equation (57) with the initial conditions [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. Table 6 analyzes the accuracy of the numerical solutions. Note: we made [Forumla omitted. See PDF.] to avoid numerical divergences at [Forumla omitted. See PDF.].

Figure 3. Solution to Equation (60) using the numerical solution explained in Section 3 (continuous gray lines) and using Equation (62) (black dashed lines) for [Forumla omitted. See PDF.] and different values of [Forumla omitted. See PDF.]. We have chosen in all cases units such that [Forumla omitted. See PDF.]. We also used [Forumla omitted. See PDF.] for both cases. In the numerical solutions [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], in the analytical solution [Forumla omitted. See PDF.]. The values of the minima and the mean square error (MSE) are shown in Table 7.

Figure 4. Solution to [Forumla omitted. See PDF.] for different values of [Forumla omitted. See PDF.] in units in which [Forumla omitted. See PDF.]. We used [Forumla omitted. See PDF.].

Figure 5. Solution to [Forumla omitted. See PDF.] for different values of [Forumla omitted. See PDF.] in units in which [Forumla omitted. See PDF.]. The continuous line shows the solution for [Forumla omitted. See PDF.], the dotted line for [Forumla omitted. See PDF.]. We used [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 6. Solution to [Forumla omitted. See PDF.] with DBC, as given in example 3 for different values. We set [Forumla omitted. See PDF.]. We made [Forumla omitted. See PDF.] in the left plot and [Forumla omitted. See PDF.] in the right one.

Figure 7. Solution to [Forumla omitted. See PDF.] with NBC, as given in example 3 for different values. We set [Forumla omitted. See PDF.] and made [Forumla omitted. See PDF.].

Table 1

A summary on the change of notation from continuous to discrete form. We have $r^0=R_{\mathrm{int}}$ and $j\in \{0,1,2,\cdots ,N\}$.

Variable or Function	Equivalent Array
r and $r^{\prime }$	$r^j=r^0+hj,$
$G_{\lambda }(r,r^{\prime },{\theta }^{\prime })$	$G_{\lambda |{\theta }^{\prime }}^{j,k}=G_{\lambda }(r^j,{r^{\prime }}^k,{\theta }^{\prime })$
$P\left(r\right)\equiv r^2$	$P^j=P\left(r^j\right)$
$Q_{\mu }\left(r\right)\equiv r^2{f}_{r\mu }\left(r\right)+r{\delta }_{0,\mu }$	$Q_{\mu }^j=Q_{\mu }\left(r^j\right)$
$R_{\lambda ,\mu }\left(r\right)\equiv r^2{g}_{\mu }\left(r\right)-{\lambda }^2{\delta }_{0,\mu }+ir(\lambda -\mu )f_{\theta \mu }\left(r\right)$	$R_{\lambda ,\mu }^j=R_{\lambda ,\mu }\left(r^j\right)$

Table 2

Non-vanishing matrix elements of $\mathbb{U}$ for $1\le j\le N-1$ for an annulus and a disc.

Mode	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }}^{\mathit{j},\mathit{j}}$	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }}^{\mathit{j},\mathit{j}\pm 1}$	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }-\mathit{\mu }}^{\mathit{j},\mathit{j}}$	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }-\mathit{\mu }}^{\mathit{j},\mathit{j}\pm 1}$
$\forall \lambda$	$2P^j-h^2{R}_{\lambda ,0}^j$	$-P^j\mp \frac{h}{2}{Q}_0^j$	$-h^2{R}_{\lambda ,\mu }^j$	$\mp \frac{h}{2}{Q}_{\mu }^j$

Table 3

Non-vanishing matrix elements of matrix $\mathbb{U}$ that define DBC (D) and NBC (N) for an annulus. N/A specifies those elements that lay outside $\mathbb{U}$.

	Mode	j	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }}^{\mathit{j},\mathit{j}}$	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }}^{\mathit{j},\mathit{j}+1}$	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }}^{\mathit{j},\mathit{j}-1}$	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }-\mathit{\mu }}^{\mathit{j},\mathit{j}}$
(D)	$\forall \lambda$	$0,N$	1	0, N/A	N/A, 0	0
(N)	$\forall \lambda$	0	$2P^0-h^2{R}_{\lambda ,0}^0$	$-2P^1$	N/A	$-h^2{R}_{\lambda ,\mu }^0$
(N)	$\forall \lambda$	N	$2P^N-h^2{R}_{\lambda ,0}^N$	N/A	$-2P^{N-1}$	$-h^2{R}_{\lambda ,\mu }^N$

Table 4

Non-vanishing matrix elements of matrix $\mathbb{U}$ that define DBC (shown above) and NBC (shown below) for a disc. N/A specifies those elements that lay outside $\mathbb{U}$.

Mode	j	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }}^{\mathit{j},\mathit{j}}$	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }}^{\mathit{j},\mathit{j}+1}$	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }}^{\mathit{j},\mathit{j}-1}$	${\mathbb{U}}_{\mathit{\lambda },\mathit{\lambda }-\mathit{\mu }}^{\mathit{j},\mathit{j}}$
$\lambda =0$	0	$2-h^2{g}_0^0$	$-2$	N/A	0
$\lambda =0$	N	1	N/A	0	0
$\lambda \ne 0$	$0,N$	1	0, N/A	N/A, 0	0
$\lambda =0$	0	$2-h^2{g}_0^0$	$-2$	N/A	0
$\forall \lambda$	N	$2P^N-h^2{R}_{\lambda ,0}^N$	N/A	$-2P^{N-1}$	$-h^2{R}_{\lambda ,\mu }^N$
$\lambda \ne 0$	0	1	0	N/A	0

Table 5

Elements ${\alpha }_{\lambda }^j$ that define DBC and NBC for an annulus (A) and a disc (D).

Mode	j	${\mathit{\alpha }}_{\mathit{\lambda }}^{\mathit{j}}$ (DBC)	${\mathit{\alpha }}_{\mathit{\lambda }}^{\mathit{j}}$ (NBC)	Geom.
$\lambda =0$	0	${ϵ}^{-1}$	${ϵ}^{-1}$	(D)
$\lambda =0$	0	0	$r^0$	(A)
$\lambda =0$	$1\le j\le N-1$	$r^j$	$r^j$	(A, D)
$\lambda =0$	N	0	$r^N$	(A, D)
$\lambda \ne 0$	0	0	0	(D)
$\lambda \ne 0$	0	0	$r^0-\frac{h^2}{4r^0}$	(A)
$\lambda \ne 0$	$1\le j\le N-1$	$r^j-\frac{h^2}{4r^j}$	$r^j-\frac{h^2}{4r^j}$	(A, D)
$\lambda \ne 0$	N	0	$r^N-\frac{h^2}{4r^N}$	(A, D)

Table 6

Analysis of the accuracy of the two numerical solutions to Equation (57). The first two blocks use the conditions $\psi \left(0\right)=0$ and $\psi \left(4\right)=2$ and show the value of the maximum—exact value ${\psi }_{max}^{\mathrm{DBC}}=2.6524822\cdots$, its percentage error (PE) (from 0 t 100%), and mean square error (MSE) of the function along the domain $[0,4]$; the first block uses the weight function, the second one does not. The block below shows something similar for the conditions ${\psi }^{\prime }\left(0\right)=0$ and ${\psi }^{\prime }\left(4\right)=4$ and focus on the global minimum—exact value ${\psi }_{min}^{\mathrm{NBC}}=-3.88574186\cdots$ using the weight function—remember that for NBC the formalism with the weight function is required.

Func.	${\mathit{\psi }}_{\mathit{N}=32}$	${\mathit{\psi }}_{\mathit{N}=128}$	${\mathit{\psi }}_{\mathit{N}=512}$
${\psi }_{max}^{\mathrm{DBC}}$	$2.6510$	$2.6526$	$2.6525$
$\mathrm{PE}$	$5.6550\times {10}^{-2}$	$3.7700\times {10}^{-3}$	$1.9221\times {10}^{-4}$
$\mathrm{MSE}$	$3.2824$	$2.6156\times {10}^{-9}$	$1.0287\times {10}^{-11}$
${\psi }_{max}^{\mathrm{DBC}}$	$2.7372$	$2.6736$	$2.6577$
$\mathrm{PE}$	$3.1939$	$7.9615\times {10}^{-1}$	$1.9671\times {10}^{-1}$
$\mathrm{MSE}$	$3.4935$	$2.1177\times {10}^{-4}$	$1.3148\times {10}^{-5}$
${\psi }_{min}^{\mathrm{NBC}}$	$-3.8852$	$-3.8853$	$-3.8856$
$\mathrm{PE}$	$1.4438\times {10}^{-2}$	$1.2060\times {10}^{-2}$	$3.1806\times {10}^{-3}$
$\mathrm{MSE}$	$1.6316\times {10}^{-3}$	$2.3096\times {10}^{-5}$	$3.5993\times {10}^{-7}$

Table 7

Values of the minima shown in Figure 3 for the numerical solution (NS) and analytical solution (AS). PE means percentage error (from 0 to $100\%$) and MSE is the mean square error.

	Values of the Minima
$G(r,r^{\prime })$	$r^{\prime }=0$	$r^{\prime }=1.172$	$r^{\prime }=3.516$	$r^{\prime }=7.031$
NS	$-1.288$	$-0.783$	$-0.609$	$-0.498$
AS	$-1.2778$	$-0.7835$	$-0.6087$	$-0.4984$
$\mathrm{PE}$	$0.8263$	$3.566\times 10-2$	$5.148\times {10}^{-3}$	$1.603\times {10}^{-3}$
$\mathrm{MSE}$	$3.629\times {10}^{-8}$	$1.343\times {10}^{-10}$	$5.053\times {10}^{-12}$	$6.379\times {10}^{-13}$

Appendix A. Deduction of the Weight Function

The relation obeyed by the weight function that makes the Liouville operator self-adjoint can be deduced by performing a direct substitution of Equation (1) into (5) and using Green’s: ${\int }_{\mathcal{V}}\varphi \left({\nabla }^2\psi \right)\mathrm{d}\mathbf{r}={\int }_{\mathcal{V}}\psi \left({\nabla }^2\varphi \right)\mathrm{d}\mathbf{r}+{\oint }_{\partial \mathcal{V}}\left[\varphi \left(\overrightarrow{\nabla }\psi \right)-\psi \left(\overrightarrow{\nabla }\varphi \right)\right]·\mathbf{n}\mathrm{d}S$ and the Divergence: ${\int }_{\mathcal{V}}\mathbf{a}·\left(\overrightarrow{\nabla }\psi \right)\mathrm{d}\mathbf{r}={\oint }_{\partial \mathcal{V}}\psi \mathbf{a}·\mathbf{n}dS-{\int }_{\mathcal{V}}\psi (\overrightarrow{\nabla }·\mathbf{a})\mathrm{d}\mathbf{r}$ theorems. After writing it conveniently, the result of this operation is (A1)$$\begin{array}{cc}\hfill \psi \left(\mathbf{r}\right)=& \hfill {\int }_{\mathbf{r}}^{\prime }G({\mathbf{r}}^{\prime },\mathbf{r})[w({\mathbf{r}}^{\prime },\mathbf{r})\left[{\nabla }_{\left\{{\mathbf{r}}^{\prime }\right\}}^2\psi \left({\mathbf{r}}^{\prime }\right)\right]\hfill \\ \hfill & \hfill +\left[{\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}w(\mathbf{r},{\mathbf{r}}^{\prime })\right]·\left[{\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}\psi \left({\mathbf{r}}^{\prime }\right)\right]+w({\mathbf{r}}^{\prime },\mathbf{r})g\left({\mathbf{r}}^{\prime }\right)\psi \left({\mathbf{r}}^{\prime }\right)]\mathrm{d}{\mathbf{r}}^{\prime }+\hfill \\ \hfill & \hfill {\int }_{{\mathbf{r}}^{\prime }}\{\left[{\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}w({\mathbf{r}}^{\prime },\mathbf{r})-w({\mathbf{r}}^{\prime },\mathbf{r})\overrightarrow{f}\left({\mathbf{r}}^{\prime }\right)\right]·{\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}\psi \left({\mathbf{r}}^{\prime }\right)\hfill \\ \hfill & \hfill +{\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}·\left[{\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}w({\mathbf{r}}^{\prime },\mathbf{r})-w({\mathbf{r}}^{\prime },\mathbf{r})\overrightarrow{f}\left({\mathbf{r}}^{\prime }\right)\right]\psi \left({\mathbf{r}}^{\prime }\right)\}\mathrm{d}{\mathbf{r}}^{\prime }+\hfill \\ \hfill & \hfill {\oint }_{\partial {\mathbf{r}}^{\prime }}[w({\mathbf{r}}^{\prime },\mathbf{r})\left[\psi \left({\mathbf{r}}^{\prime }\right){\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}G({\mathbf{r}}^{\prime },\mathbf{r})-G({\mathbf{r}}^{\prime },\mathbf{r}){\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}\psi \left({\mathbf{r}}^{\prime }\right)\right]\hfill \\ \hfill & \hfill -G({\mathbf{r}}^{\prime },\mathbf{r})\psi \left({\mathbf{r}}^{\prime }\right)\left[{\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}w({\mathbf{r}}^{\prime },\mathbf{r})-w({\mathbf{r}}^{\prime },\mathbf{r})\overrightarrow{f}\left({\mathbf{r}}^{\prime }\right)\right]]·\mathbf{n}d{S}^{\prime }.\hfill \end{array}$$

Notice that under the choice ${\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}w({\mathbf{r}}^{\prime },\mathbf{r})-w({\mathbf{r}}^{\prime },\mathbf{r})\overrightarrow{f}\left({\mathbf{r}}^{\prime }\right)=0$, last equation transforms into Equation (8).

Appendix B. Finite Elements Method, Matrix Elements

The elements of matrices ${\mathbf{A}}^{\left(n\right)}$ introduced in Section 3.1 depend on the required level of accuracy and the central site $\eta$ that we choose; a general algorithm to deduce such elements is shown in [42]. In the simplest case, as we choose $\eta =0$ as the central point in conjunction with the two closets neighbors—three-point stencil approximation—we have the following relations [42,43] (A2)$$\begin{array}{cc}\hfill {\mathbf{A}}^{\left(1\right)}& \hfill ={\left(\begin{array}{ccccccc}\hfill 0& \hfill 1& \hfill 0& \hfill \cdots & \hfill 0& \hfill 0& \hfill 0\\ \hfill -1& \hfill 0& \hfill 1& \hfill \cdots & \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill -1& \hfill 0& \hfill \cdots & \hfill 0& \hfill 0& \hfill 0\\ \hfill ⋮& \hfill ⋮& \hfill ⋮& \hfill \ddots & \hfill ⋮& \hfill ⋮& \hfill ⋮\\ \hfill 0& \hfill 0& \hfill 0& \hfill \cdots & \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill \cdots & \hfill -1& \hfill 0& \hfill 1\\ \hfill 0& \hfill 0& \hfill 0& \hfill \cdots & \hfill 0& \hfill -1& \hfill 0\end{array}\right)}_{N+1\times N+1},\hfill \end{array}$$ (A3)$$\begin{array}{cc}\hfill {\mathbf{A}}^{\left(2\right)}& \hfill ={\left(\begin{array}{ccccccc}\hfill -2& \hfill 1& \hfill 0& \hfill \cdots & \hfill 0& \hfill 0& \hfill 0\\ \hfill 1& \hfill -2& \hfill 1& \hfill \cdots & \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill -2& \hfill \cdots & \hfill 0& \hfill 0& \hfill 0\\ \hfill ⋮& \hfill ⋮& \hfill ⋮& \hfill \ddots & \hfill ⋮& \hfill ⋮& \hfill ⋮\\ \hfill 0& \hfill 0& \hfill 0& \hfill \cdots & \hfill -2& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill \cdots & \hfill 1& \hfill -2& \hfill 1\\ \hfill 0& \hfill 0& \hfill 0& \hfill \cdots & \hfill 0& \hfill 1& \hfill -2\end{array}\right)}_{N+1\times N+1}.\hfill \end{array}$$

Notice how the matrices ${\mathbf{A}}^{\left(n\right)}$ must be truncated at the boundaries; this is a natural consequence of the FEM, coming from the boundary conditions.

Appendix C. Derivatives at the Boundaries for DBC

If the Green function is used to find a non-homgeneous function with DBC, the derivatives of the Green function at the borders are needed—see Equations (8) and (14). Combining Equation (44) with the boundary conditions stated in Equations (37), (39), (41) and (43), we deduce the following two relations ($j=0$ and $j=N$ refer to the two possible boundaries) (A4)$$\begin{array}{cc}\hfill {G^{\prime }}_{\theta ,{\theta }^{\prime }}^{(0,N),k}& \hfill =\frac{1}{2\pi }\sum _{\lambda ,\mu }{e}^{i\lambda \theta }{e}^{-i\mu {\theta }^{\prime }}{\mathbb{B}}_{\lambda ,\mu }^{(0,N),k},\mathrm{where}\hfill \\ \hfill {\mathbb{B}}_{\lambda ,\mu }^{(0,N),k}& \hfill =\mp P^{0,N}\sum _{\nu }{\left[{\mathbb{S}}^{(0,N)}\right]}_{\lambda ,\nu }^{-1}{\mathbb{A}}_{\nu ,\mu }^{(1,N-1),k}.\hfill \end{array}$$

The matrix elements associated with ${\mathbb{S}}^{(0,N)}$ are ${\mathbb{S}}_{\lambda ,\lambda }^{(0,N)}=P^{0,N}\mp \frac{h}{2}{Q}_0^{0,N}$ and ${\mathbb{S}}_{\lambda ,\lambda -\mu }^{(0,N)}=\mp \frac{h}{2}{Q}_{\mu }^{0,N}$. Since $\psi$ is known at the boundaries, the derivatives ${G^{\prime }}_{\theta ,{\theta }^{\prime }}^{0,0}$ and ${G^{\prime }}_{\theta ,{\theta }^{\prime }}^{N,N}$ are irrelevant; additionally, a disk only requires the calculation of ${G^{\prime }}_{\theta ,{\theta }^{\prime }}^{N,k}$. The symmetric elements ${G^{\prime }}_{\theta ,{\theta }^{\prime }}^{j(0,N)}$ can be found similarly, in terms of the transpose elements ${\mathbb{B}}_{\lambda ,\mu }^{j,(0,N)}$, which are defined according to last expression by performing the index change and transposition ${\mathbb{A}}_{\nu ,\mu }^{(1,N-1),k}\to {\mathbb{A}}_{\nu ,\mu }^{j,(1,N-1)}$. In one dimension the derivatives are ${G^{\prime }}^{(0,N)k}=\mp \frac{P^{0,N}{\mathbb{A}}^{(1,N-1),k}}{P^{0,N}\mp \frac{h}{2}{Q}^{0,N}}$.

Appendix D. Symmetry Properties of Some Matrix Elements

Expanding Equation (48) to eliminate the negative modes, we can write the Green function as (A5)$$\begin{array}{cc}\hfill G_{\theta ,{\theta }^{\prime }}^{j,k}& =-\frac{h}{2\pi }[{\mathbb{A}}_{0,0}^{j,k}+\sum _{\lambda \ge 1}\left\{\left[{\mathbb{A}}_{-\lambda ,0}^{j,k}+{\mathbb{A}}_{\lambda ,0}^{j,k}\right]cos\left(\lambda \theta \right)+i\left[-{\mathbb{A}}_{-\lambda ,0}^{j,k}+{\mathbb{A}}_{\lambda ,0}^{j,k}\right]sin\left(\lambda \theta \right)\right\}\hfill \\ \hfill & \phantom{=}+\sum _{\mu \ge 1}\left\{\left[{\mathbb{A}}_{0,-\mu }^{j,k}+{\mathbb{A}}_{0,\mu }^{j,k}\right]cos\left(\mu {\theta }^{\prime }\right)+i\left[{\mathbb{A}}_{0,-\mu }^{j,k}-{\mathbb{A}}_{0,\mu }^{j,k}\right]sin\left(\mu \theta \right)\right\}\hfill \\ \hfill & \phantom{=}+\sum _{\lambda ,\mu \ge 1}\{\left[{\mathbb{A}}_{-\lambda ,-\mu }^{j,k}+{\mathbb{A}}_{-\lambda ,\mu }^{j,k}+{\mathbb{A}}_{\lambda ,-\mu }^{j,k}+{\mathbb{A}}_{\lambda ,\mu }^{j,k}\right]cos\left(\lambda \theta \right)cos\left(\mu {\theta }^{\prime }\right)\hfill \\ \hfill & \phantom{=}+i\left[{\mathbb{A}}_{\lambda ,-\mu }^{j,k}+{\mathbb{A}}_{\lambda ,\mu }^{j,k}-{\mathbb{A}}_{-\lambda ,-\mu }^{j,k}-{\mathbb{A}}_{-\lambda ,\mu }^{j,k}\right]sin\left(\lambda \theta \right)cos\left(\mu {\theta }^{\prime }\right)\hfill \\ \hfill & \phantom{=}+i\left[{\mathbb{A}}_{-\lambda ,-\mu }^{j,k}-{\mathbb{A}}_{-\lambda ,\mu }^{j,k}+{\mathbb{A}}_{\lambda ,-\mu }^{j,k}-{\mathbb{A}}_{\lambda ,\mu }^{j,k}\right]cos\left(\lambda \theta \right)sin\left(\mu {\theta }^{\prime }\right)\hfill \\ \hfill & \phantom{=}+\left[{\mathbb{A}}_{-\lambda ,-\mu }^{j,k}-{\mathbb{A}}_{-\lambda ,\mu }^{j,k}-{\mathbb{A}}_{\lambda ,-\mu }^{j,k}+{\mathbb{A}}_{\lambda ,\mu }^{j,k}\right]sin\left(\lambda \theta \right)sin\left(\mu {\theta }^{\prime }\right)\left\}\right].\hfill \end{array}$$

Since the Green function must be real for real Liouville operators, we demand that the imaginary contributions of last expression must vanish. Hence, we have the restrictions stated in Equation (49).

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 Green function is real. Using the properties stated in Equation (49) into Equation (A5), we find that (A6)$$\begin{array}{cc}\hfill G_{\theta {\theta }^{\prime }}^{jk}=& \hfill -\frac{h}{2\pi }\mathrm{Re}\left({\mathbb{A}}_{0,0}^{j,k}\right)-\frac{h}{\pi }\sum _{\lambda \ge 1}[\mathrm{Re}\left({\mathbb{A}}_{\lambda ,0}^{j,k}\right)cos\left(\lambda \theta \right)+\mathrm{Re}\left({\mathbb{A}}_{0,\lambda }^{j,k}\right)cos\left(\lambda {\theta }^{\prime }\right)\hfill \\ \hfill & \hfill -\mathrm{Im}\left({\mathbb{A}}_{\lambda ,0}^{j,k}\right)sin\left(\lambda \theta \right)+\mathrm{Im}\left({\mathbb{A}}_{0,\lambda }^{j,k}\right)sin\left(\lambda {\theta }^{\prime }\right)]\hfill \\ \hfill & \hfill -\frac{h}{\pi }\sum _{\lambda ,\mu \ge 1}\left[\mathrm{Re}\left({\mathbb{A}}_{\lambda ,\mu }^{j,k}\right)cos(\lambda \theta -\mu {\theta }^{\prime })+\mathrm{Re}\left({\mathbb{A}}_{\lambda ,-\mu }^{j,k}\right)cos(\lambda \theta +\mu {\theta }^{\prime })\right]\hfill \\ \hfill & \hfill +\frac{h}{\pi }\sum _{\lambda ,\mu \ge 1}\left[\mathrm{Im}\left({\mathbb{A}}_{\lambda ,\mu }^{j,k}\right)sin(\lambda \theta -\mu {\theta }^{\prime })+\mathrm{Im}\left({\mathbb{A}}_{\lambda ,-\mu }^{j,k}\right)sin(\lambda \theta +\mu {\theta }^{\prime })\right].\hfill \end{array}$$

In the presence of angular symmetry, $A_{\lambda \mu }^{jk}=A_{\lambda \mu }^{jk}{\delta }_{\lambda \mu }$, so last equation reduces to (A7)$$\begin{array}{cc}\hfill G_{\theta {\theta }^{\prime }}^{jk}=& -\frac{h}{2\pi }\mathrm{Re}\left({\mathbb{A}}_{0,0}^{j,k}\right)\hfill \\ \hfill & -\frac{h}{\pi }\sum _{\lambda \ge 1}\left[\mathrm{Re}\left({\mathbb{A}}_{\lambda ,\lambda }^{j,k}\right)cos\left[\lambda (\theta -{\theta }^{\prime })\right]-\mathrm{Im}\left({\mathbb{A}}_{\lambda ,\lambda }^{j,k}\right)sin\left[\lambda (\theta -{\theta }^{\prime })\right]\right].\hfill \end{array}$$

Appendix F. Computation of ψ(r) as an Exponential Expansion

In this section, we will derive expressions for Equations (8) and (14) for both DBC and NBC. Although both approaches must lead to the same results, it is worthwhile to show how both relations can be found through the formalism we have described.

For convenience, we will split $\psi (r,\theta )$ into a volume (V) and surface (S) contribution—the volume contribution is the term containing the integral over ${\mathbf{r}}^{\prime }$ in Equations (8) and (14); the surface contribution is the one containing the integral over the closed surface $\partial {\mathbf{r}}^{\prime }$ in the same equations. For DBC and NBC, $\psi$ can be written in discrete coordinates as (A8a)$$\begin{array}{cc}\hfill {\left({\psi }_{\theta }^j\right)}^{\mathrm{DBC}}& \hfill ={\left({\psi }_{\theta }^j\right)}_V+{\left({\psi }_{\theta }^j\right)}_S^{\mathrm{DBC}},\hfill \end{array}$$ (A8b)$$\begin{array}{cc}\hfill {\left({\psi }_{\theta }^j\right)}^{\mathrm{NBC}}& \hfill ={\left({\psi }_{\theta }^j\right)}_V+{\left({\psi }_{\theta }^j\right)}_S^{\mathrm{NBC}}.\hfill \end{array}$$

There are many ways to evaluate numerically an integral. We will use one of the simplest, however, very efficient, ways to do so, the so called trapezoid rule. Due to the discretization we have used, this rule will be applied to evaluate the the radial integrals, appearing in the volume contributions; the integrals over angular coordinates will be evaluated directly using the Fourier expansions of the functions involved.

• Using the weight function

Having adopted the convention described in Equation (8), we start performing a Fourier expansions of the external field: $\varphi \left({\mathbf{r}}^{\prime }\right)\to {\varphi }_{{\theta }^{\prime }}^k={\sum }_{\lambda }{\varphi }_{\lambda }^k{e}^{i\lambda {\theta }^{\prime }}$, the weight function: $w({\mathbf{r}}^{\prime },\mathbf{r})\to \frac{w_{{\theta }^{\prime }}^k}{w(r^j,\theta )}=\frac{1}{w(r^j,\theta )}{\sum }_{\lambda }{w}_{\lambda }^k{e}^{i\lambda {\theta }^{\prime }}$—and something similar for the boundary conditions ${\psi }_{{\theta }^{\prime }}^{(0,N)}$ and ${{\psi }^{\prime }}_{{\theta }^{\prime }}^{(0,N)}$. Now, we will define the function (A9)$$\begin{array}{cc}\hfill {\xi }_w(k,M^{k,j},{\eta }^k,\theta )& \hfill =\frac{1}{2\pi }{\int }_0^{2\pi }d{\theta }^{\prime }\sum _{\lambda ,\mu }{e}^{i\lambda {\theta }^{\prime }}{M}_{\lambda ,\mu }^{k,j}{e}^{-i\mu \theta }\sum _{\nu }{\eta }_{\nu }^k{e}^{i\nu {\theta }^{\prime }}\sum _{\rho }{w}_{\rho }^k{e}^{i\rho {\theta }^{\prime }}\hfill \\ \hfill & \hfill =\sum _{\lambda ,\mu ,\nu }{e}^{-i\mu \theta }{M}_{\lambda ,\mu }^{k,j}{\eta }_{\nu }^k{w}_{-\lambda -\nu }^k.\hfill \end{array}$$

This definition will be used to define the volume- and surface-terms.

Since the volume-term can be written as ${\int }_{r^0}^{r^N}{r}^{\prime }d{r}^{\prime }{\int }_0^{2\pi }w({\mathbf{r}}^{\prime },\mathbf{r})G({\mathbf{r}}^{\prime },\mathbf{r})\varphi \left({\mathbf{r}}^{\prime }\right)d{\theta }^{\prime }$ ($r^0=R_{\mathrm{int}},r^N=R_{\mathrm{ext}}$), we can say that (A10)$$\begin{array}{cc}\hfill {\left({\psi }_{\theta }^j\right)}_V=& \hfill -\frac{h^2}{w(r^j,\theta )}[\sum _{k=1}^{N-1}{r}^k{\xi }_w(k,{\mathbb{A}}^{k,j},{\varphi }^k,\theta )+\frac{1}{2}{r}^0{\xi }_w(0,{\mathbb{A}}^{0,j},{\varphi }^0,\theta )\hfill \\ \hfill & \hfill +\frac{1}{2}{r}^N{\xi }_w(N,{\mathbb{A}}^{N,j},{\varphi }^0,\theta )].\hfill \end{array}$$

The surface-term arising in DBC is written as $r^{\prime }{\int }_0^{2\pi }w({\mathbf{r}}^{\prime },\mathbf{r})\psi \left({\mathbf{r}}^{\prime }\right){\partial }_{r^{\prime }}G({\mathbf{r}}^{\prime },\mathbf{r})d{\theta }^{\prime }{|}_{r^0}^{r^N}$. Similarly as shown above, in discrete coordinates it is given by (A11)$$\begin{array}{cc}\hfill {\left({\psi }_{\theta }^j\right)}_S^{\mathrm{DBC}}& \hfill =\frac{1}{w(r^j,\theta )}\left[r^N{\xi }_w(N,{\mathbb{B}}^{N,j},{\psi }^N,\theta )-r^0{\xi }_w(0,{\mathbb{B}}^{0,j},{\psi }^0,\theta )\right].\hfill \end{array}$$

Finally, the surface-term $-{\oint }_{\partial {\mathbf{r}}^{\prime }}w({\mathbf{r}}^{\prime },\mathbf{r})G({\mathbf{r}}^{\prime },\mathbf{r}){\overrightarrow{\nabla }}_{\left\{{\mathbf{r}}^{\prime }\right\}}\psi \left({\mathbf{r}}^{\prime }\right)·\mathbf{n}d{S}^{\prime }$ that appears in NBC can now be expanded in the form $-r^{\prime }{\int }_0^{2\pi }w({\mathbf{r}}^{\prime },\mathbf{r})G({\mathbf{r}}^{\prime },\mathbf{r}){\partial }_{r^{\prime }}\psi \left({\mathbf{r}}^{\prime }\right)d{\theta }^{\prime }{|}_{r^0}^{r^N}$; it now becomes (A12)$$\begin{array}{cc}\hfill {\left({\psi }_{\theta }^j\right)}_S^{\mathrm{NBC}}& \hfill =\frac{h}{w(r^j,\theta )}\left[r^N{\xi }_w(N,{\mathbb{A}}^{N,j},{\psi }^N,\theta )-r^0{\xi }_w(0,{\mathbb{A}}^{0,j},{\psi }^0,\theta )\right].\hfill \end{array}$$

Remark A1.

The matrix elements of matrices $\mathbb{A}$ and ${\mathbb{B}}^{(0,N)}$ are given by Equations (48) and Equation (A4), respectively—the indices associated to the position in the blocks have been omitted by convenience. The function ${\psi }_{\theta }^j$ is defined in the interval $1\le j\le N-1$; in DBC the terms ${\psi }_{\theta }^0$ and ${\psi }_{\theta }^N$ are given, in NBC the function at the borders is not accurate enough due to the discontinuity of the Green function at the borders.

• Using no weight function

When we adapt the convention stated in Equation (14), Equations (A10)–(A12) are slightly modified. We now define the function $\xi$ as (A13)$$\begin{array}{cc}\hfill \xi (k,M^{j,k},{\eta }^k,\theta )=& \hfill {\int }_0^{2\pi }\frac{d{\theta }^{\prime }}{2\pi }\sum _{\lambda ,\mu }{e}^{i\lambda \theta }{M}_{\lambda ,\mu }^{j,k}{e}^{-i\mu {\theta }^{\prime }}\sum _{\nu }{\eta }_{\nu }^k{e}^{i\nu {\theta }^{\prime }}=\sum _{\lambda ,\mu }{e}^{i\lambda \theta }{M}_{\lambda ,\mu }^{j,k}{\eta }_{\mu }^k.\hfill \end{array}$$

We can now conclude that (A14)$$\begin{array}{cc}\hfill {\left({\psi }_{\theta }^j\right)}_V& \hfill =-h^2\sum _{k=1}^{N-1}[r^k\xi (k,{\mathbb{A}}^{j,k},{\varphi }^k,\theta )+\frac{1}{2}{r}^0\xi (0,{\mathbb{A}}^{j,0},{\varphi }^0,\theta )\hfill \end{array}$$ (A15)$$\begin{array}{cc}\hfill & \hfill +\frac{1}{2}{r}^N\xi (N,{\mathbb{A}}^{j,N},{\varphi }^N,\theta )],\hfill \end{array}$$ (A16)$$\begin{array}{cc}\hfill {\left({\psi }_{\theta }^j\right)}_S^{\mathrm{DBC}}& \hfill =r^N\xi (N,{\mathbb{B}}^{j,N},{\psi }^N,\theta )-r^0\xi (0,{\mathbb{B}}^{j,0},{\psi }^0,\theta ),\hfill \\ \hfill {\left({\psi }_{\theta }^j\right)}_S^{\mathrm{NBC}}& \hfill =h{r}^N\left[\xi (N,{\mathbb{A}}^{j,N},f_r^N,\theta )+\xi (N,{\mathbb{A}}^{j,N},{{\psi }^{\prime }}^N,\theta )\right]\hfill \end{array}$$ (A17)$$\begin{array}{cc}\hfill & \hfill \phantom{=}-h{r}^0\left[\xi (0,{\mathbb{A}}^{j,0},f_r^0,\theta )+\xi (0,{\mathbb{A}}^{j,0},{{\psi }^{\prime }}^0,\theta )\right].\hfill \end{array}$$

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 the expressions found are much longer, this allows us to use the symmetry properties, Equation (49), to get rid of irrelevant terms and explicitly express $\psi \left(\mathbf{r}\right)$ as a real function. In addition, the exponential expansion defined in Appendix F might introduce some spurious imaginary contributions, which may arise by as a consequence of the truncating process of matrix $\mathbb{U}$—the complex conjugate counterparts of some modes may be discarded in this process. Taking advantage of the definitions used in Appendix F, Equations (A10) to (A12) are still valid when we adopt the convention stated in Equation (8); similarly, when the convention Equation (14) is adopted, Equations (A14) to (A17) are also valid. Now, we only need to expand ${\xi }_w$ and $\xi$ eliminating the negative complex modes to express them as sum of real modes. By doing so, Equation (A9) becomes (A18)$$\begin{array}{cc}\hfill & {\xi }_w(k,M^{k,j},{\eta }^k,\theta )=\mathrm{Re}\left(M_{0,0}^{k,j}\right)\mathrm{Re}\left({\eta }_0^k\right)\mathrm{Re}\left(w_0^k\right)+2\sum _{\lambda \ge 1}\mathrm{Re}\left({\eta }_0^k\right)\left[\mathrm{Re}\left(M_{0,0}^{k,j}\right)X_{\lambda }^k+Y_{\lambda }^{k,j}\right]\hfill \\ \hfill & +2\sum _{\lambda \ge 1}\mathrm{Re}\left({\eta }_0^k\right)\mathrm{Re}\left(w_0^k\right)\left[\mathrm{Re}\left(M_{0,\lambda }^{k,j}\right)cos\left(\lambda \theta \right)+\mathrm{Im}\left(M_{0,\lambda }^{k,j}\right)sin\left(\lambda \theta \right)\right]\hfill \\ \hfill & +4\sum _{\lambda ,\mu \ge 1}{X}_{\mu }^k\left[\mathrm{Re}\left(M_{0,\lambda }^{k,j}\right)cos\left(\lambda \theta \right)+\mathrm{Im}\left(M_{0,\lambda }^{k,j}\right)sin\left(\lambda \theta \right)\right]\hfill \\ \hfill & +2\sum _{\lambda ,\mu \ge 1}\mathrm{Re}\left({\eta }_0^k\right)\left[\mathrm{Re}\left(w_{\lambda }^k\right)M_{(+)\lambda ,\mu }^{\left(1\right)k,j}+\mathrm{Im}\left(w_{\lambda }^k\right)M_{(+)\lambda ,\mu }^{\left(2\right)k,j}\right]cos\left(\mu \theta \right)\hfill \\ \hfill & -2\sum _{\lambda ,\mu \ge 1}\mathrm{Re}\left({\eta }_0^k\right)\left[\mathrm{Im}\left(w_{\lambda }^k\right)M_{(-)\lambda ,\mu }^{\left(1\right)k,j}-\mathrm{Re}\left(w_{\lambda }^k\right)M_{(-)\lambda ,\mu }^{\left(2\right)k,j}\right]sin\left(\mu \theta \right)\hfill \\ \hfill & +2\sum _{\lambda ,\mu \ge 1}\mathrm{Re}\left({\eta }_{\mu }^k\right)\left[\mathrm{Re}\left(M_{\lambda ,0}^{k,j}\right)w_{+(\lambda ,\mu )}^{\left(1\right)k}+\mathrm{Im}\left(M_{\lambda ,0}^{k,j}\right)w_{+(\lambda ,\mu )}^{\left(2\right)k}\right]\hfill \\ \hfill & +2\sum _{\lambda ,\mu \ge 1}\mathrm{Im}\left({\eta }_{\mu }^k\right)\left[\mathrm{Re}\left(M_{\lambda ,0}^{k,j}\right)w_{-(\lambda ,\mu )}^{\left(2\right)k}-\mathrm{Im}\left(M_{\lambda ,0}^{k,j}\right)w_{-(\lambda ,\mu )}^{\left(1\right)k}\right]\hfill \\ \hfill & +2\sum _{\lambda ,\mu ,\nu \ge 1}\left[M_{(+)\lambda ,\mu }^{\left(1\right)k,j}\left[Z_{(\lambda ,\nu )}^{\left(1\right)k}+Z_{(\lambda ,\nu )}^{\left(2\right)k}\right]+M_{(+)\lambda ,\mu }^{\left(2\right)k,j}\left[Z_{(\lambda ,\nu )}^{\left(3\right)k}+Z_{(\lambda ,\nu )}^{\left(4\right)k}\right]\right]cos\left(\mu \theta \right)\hfill \\ \hfill & +2\sum _{\lambda ,\mu ,\nu \ge 1}\left[M_{(-)\lambda ,\mu }^{\left(2\right)k,j}\left[Z_{(\lambda ,\nu )}^{\left(1\right)k}+Z_{(\lambda ,\nu )}^{\left(2\right)k}\right]-M_{(-)\lambda ,\mu }^{\left(1\right)k,j}\left[Z_{(\lambda ,\nu )}^{\left(3\right)k}-Z_{(\lambda ,\nu )}^{\left(4\right)k}\right]\right]sin\left(\mu \theta \right),\hfill \end{array}$$ where we used the definitions (A19)$$\begin{array}{cc}\hfill & \hfill X_{\lambda }^k=\mathrm{Re}\left({\eta }_{\lambda }^k\right)\mathrm{Re}\left(w_{\lambda }^k\right)+\mathrm{Im}\left({\eta }_{\lambda }^k\right)\mathrm{Im}\left(w_{\lambda }^k\right)\hfill \end{array}$$ (A20)$$\begin{array}{cc}\hfill & \hfill Y_{\lambda }^{k,j}=\mathrm{Re}\left(M_{\lambda ,0}^{k,j}\right)\mathrm{Re}\left(w_{\lambda }^k\right)+\mathrm{Im}\left(M_{\lambda ,0}^{k,j}\right)\mathrm{Im}\left(w_{\lambda }^k\right);\hfill \end{array}$$ (A21)$$\begin{array}{cc}\hfill & \hfill M_{(\pm )\lambda ,\mu }^{\left(1\right)k,j}=\mathrm{Re}\left(M_{\lambda ,\mu }^{k,j}\right)\pm \mathrm{Re}\left(M_{\lambda ,-\mu }^{k,j}\right),M_{(\pm )\lambda ,\mu }^{\left(2\right)k,j}=\mathrm{Re}\left(M_{\lambda ,\mu }^{k,j}\right)\pm \mathrm{Re}\left(M_{\lambda ,-\mu }^{k,j}\right);\hfill \end{array}$$ (A22)$$\begin{array}{cc}\hfill & \hfill w_{\pm (\lambda ,\nu )}^{\left(1\right)k}=\mathrm{Re}\left(w_{\lambda +\nu }^k\right)\pm \mathrm{Re}\left(w_{\lambda -\nu }^k\right),w_{\pm (\lambda ,\nu )}^{\left(2\right)k}=\mathrm{Im}\left(w_{\lambda +\nu }^k\right)\pm \mathrm{Im}(w_{\lambda -\nu }^k);\hfill \\ \hfill & \hfill Z_{(\lambda ,\nu )}^{\left(1\right)k}=\mathrm{Re}\left({\eta }_{\nu }^k\right)w_{+(\lambda ,\nu )}^{\left(1\right)k},Z_{(\lambda ,\nu )}^{\left(2\right)k}=\mathrm{Im}\left({\eta }_{\nu }^k\right)w_{-(\lambda ,\nu )}^{\left(2\right)k},\hfill \end{array}$$ (A23)$$\begin{array}{cc}\hfill & \hfill Z_{(\lambda ,\nu )}^{\left(3\right)k}=\mathrm{Re}\left({\eta }_{\nu }^k\right)w_{+(\lambda ,\nu )}^{\left(2\right)k},Z_{(\lambda ,\nu )}^{\left(4\right)k}=\mathrm{Im}\left({\eta }_{\nu }^k\right)w_{-(\lambda ,\nu )}^{\left(1\right)k}.\hfill \end{array}$$

Similarly, Equation (A13) is now written as (A24)$$\begin{array}{cc}\hfill \xi (k,M^{j,k},{\eta }^k,\theta )& =\mathrm{Re}\left(M_{0,0}^{j,k}\right)\mathrm{Re}\left({\eta }_0^k\right)+\hfill \\ \hfill & \phantom{=}2\sum _{\lambda \ge 1}\left[R_{\lambda }^{j,k}+\mathrm{Re}\left({\eta }_0^k\right)\left[\mathrm{Re}\left(M_{\lambda ,0}^{j,k}\right)cos\left(\lambda \theta \right)-\mathrm{Im}\left(M_{\lambda ,0}^{j,k}\right)sin\left(\lambda \theta \right)\right]\right]+\hfill \\ \hfill & \phantom{=}2\sum _{\lambda ,\mu \ge 1}\left[M_{(+)\lambda ,\mu }^{\left(1\right)j,k}\mathrm{Re}\left({\eta }_{\mu }^k\right)-M_{(-)\lambda ,\mu }^{\left(2\right)j,k}\mathrm{Im}\left({\eta }_{\mu }^k\right)\right]cos\left(\lambda \theta \right)-\hfill \\ \hfill & \phantom{=}2\sum _{\lambda ,\mu \ge 1}\left[M_{(+)\lambda ,\mu }^{\left(2\right)j,k}\mathrm{Re}\left({\eta }_{\mu }^k\right)+M_{(-)\lambda ,\mu }^{\left(1\right)k,j}\mathrm{Im}\left({\eta }_{\mu }^k\right)\right]sin\left(\lambda \theta \right),\hfill \end{array}$$ with $R_{\lambda }^{j,k}=\mathrm{Re}\left(M_{0,\lambda }^{j,k}\right)\mathrm{Re}\left({\eta }_{\lambda }^k\right)-\mathrm{Im}\left(M_{0,\lambda }^{j,k}\right)\mathrm{Im}\left({\eta }_{\lambda }^{\mathbb{I}k}\right)$.

• The one-dimensional case

The function $\psi$, satisfying the equation ${\mathcal{L}}_x\psi \left(x\right)=\varphi \left(x\right)$, where ${\mathcal{L}}_x$ is defined according to Equation (52), can be found by means of the relations below. For DBC with either weight function or not (A25a)$$\begin{gathered} \begin{array}{cc}\hfill {\psi }^j& =-h^2\left[\sum _{k=1}^{N-1}{w}^{k,j}{A}^{k,j}{\varphi }^k+\frac{1}{2}{w}^{0,j}{A}^{0,j}{\varphi }^0+\frac{1}{2}{w}^{N,j}{A}^{N,j}{\varphi }^N\right]\hfill \end{array} \\ \begin{array}{cc}\hfill & \phantom{=}+\frac{P^N{w}^{N,j}{\psi }^N{A}^{N-1,j}}{1+\frac{h}{2}\frac{1}{Q^N}}+\frac{P^0{w}^{0,j}{\psi }^0{A}^{1,j}}{1-\frac{h}{2}\frac{1}{Q^0}};\hfill \end{array} \end{gathered}$$ (A25b)$$\begin{array}{cc}\hfill & \hfill =-h^2\left[\sum _{k=1}^{N-1}{A}^{j,k}{\varphi }^k+\frac{1}{2}{A}^{j,0}{\varphi }^0+\frac{1}{2}{A}^{j,N}{\varphi }^N\right]+\frac{q^N{\psi }^N{A}^{j,N-1}}{1+\frac{h}{2}\frac{1}{Q^N}}+\frac{P^0{\psi }^0{A}^{j,1}}{1-\frac{h}{2}\frac{1}{Q^0}}.\hfill \end{array}$$

For NBC (A26a)$$\begin{array}{cc}\hfill {\psi }^j& \hfill =-h^2\left[\sum _{k=1}^{N-1}{w}^{k,j}{A}^{k,j}{\varphi }^k+\frac{1}{2}{w}^{0,j}{A}^{0,j}{\varphi }^0+\frac{1}{2}{w}^{N,j}{A}^{Nj}{\varphi }^N\right]+h{P}^k{w}^{k,j}{{\psi }^{\prime }}^k{A}^{k,j}{|}_{k=0}^N;\hfill \end{array}$$ (A26b)$$\begin{array}{cc}\hfill & \hfill =-h^2\left[\sum _{k=1}^{N-1}{A}^{j,k}{\varphi }^k+\frac{1}{2}{A}^{j,0}{\varphi }^0+\frac{1}{2}{A}^{j,N}{\varphi }^N\right]+h{P}^k{A}^{j,k}\left[{{\psi }^{\prime }}^k+{\psi }^k{f}^k\right]{|}_{k=0}^N.\hfill \end{array}$$

Appendix H. Notation and List of Used Parameters

Table A1 shows a list of most parameters that have been used in this document specifying their notation

Table A1

Notation used in this document. For additional information, please see Table 1.

Parameter	Notation	Parameter	Notation
Liouville operator	${\widehat{\mathcal{L}}}_{\left\{\mathbf{r}\right\}}$	Arbitrary scalar field	□
Vector field	$\overrightarrow{f}\left(\mathbf{r}\right)$	Scalar field	$g\left(\mathbf{r}\right)$
Unknown scalar function	$\psi \left(\mathbf{r}\right)$	Non-homogeneous scalar function	$\varphi \left(\mathbf{r}\right)$
Green function	$G(\mathbf{r},{\mathbf{r}}^{\prime })$	Weight function	$w(\mathbf{r},{\mathbf{r}}^{\prime })$
Internal radius	$R_{\mathrm{int}},r^0$	External radius	$R_{\mathrm{ext}},r^N$
Polar coordinates	$r,\theta ,r^{\prime },{\theta }^{\prime }$	Fourier modes	Indexed by $\lambda$ or $\mu$
Min. radius	$ϵ$ (cutoff)	Max. Fourier mode	L (cutoff)
Grid size	h	Matrix size	N
Full “Green matrix”	${\mathbb{G}}_{{\theta }^{\prime }}$	Blocks of ${\mathbb{G}}_{{\theta }^{\prime }}$	${\mathbf{G}}_{\lambda |{\theta }^{\prime }}$
Convenient matrices	$\mathbf{P}$, ${\mathbf{Q}}_{\mu }$, ${\mathbf{R}}_{\lambda ,\mu }$	Matrix to be inverted	$\mathbb{U}$
Temporary matrices	$\mathbb{V},{\mathbb{E}}_{{\theta }^{\prime }}$, $\mathbb{A}={\mathbb{U}}^{-1}·\mathbb{V}$	Matrices used in FEM	${\mathbf{A}}^{\left(n\right)}$

References

1. Schwinger, J. Brownian motion of a quantum oscillator. J. Math. Phys.; 1961; 2, pp. 407-432. [DOI: https://dx.doi.org/10.1063/1.1703727]

2. Wang, J.S.; Agarwalla, B.K.; Li, H.; Thingna, J. Nonequilibrium green’s function method for quantum thermal transport. Front. Phys.; 2014; 9, pp. 673-697. [DOI: https://dx.doi.org/10.1007/s11467-013-0340-x]

3. Foster, S.; Neophytou, N. Effectiveness of nanoinclusions for reducing bipolar effects in thermoelectric materials. Comput. Mater. Sci.; 2019; 164, pp. 91-98. [DOI: https://dx.doi.org/10.1016/j.commatsci.2019.04.005]

4. Moulhim, A.; Tripathi, B.; Kumar, M. Nonequilibrium green function technique for analyzing electron transport through single and two levels of interacting quantum dot. Phys. Scr.; 2021; 96, 125802. [DOI: https://dx.doi.org/10.1088/1402-4896/ac198f]

5. Kadanoff, L.P.; Baym, G. Quantum Statistical Mechanics: Green’s Function Methods in Equilibrium and Nonequilibrium Problems; W. A. Benhamin, Inc. XI: New York, NY, USA, 1962; 203.

6. Hidaka, Y.; Pu, S.; Wang, Q.; Yang, D. Foundations and applications of quantum kinetic theory. Prog. Part. Nucl. Phys.; 2022; 127, 103989. [DOI: https://dx.doi.org/10.1016/j.ppnp.2022.103989]

7. Alkofer, R.; von Smekal, L. The infrared behaviour of qcd green’s functions: Confinement, dynamical symmetry breaking, and hadrons as relativistic bound states. Phys. Rep.; 2001; 353, pp. 281-465. [DOI: https://dx.doi.org/10.1016/S0370-1573(01)00010-2]

8. Lucarini, V. Revising and extending the linear response theory for statistical mechanical systems: Evaluating observables as predictors and predictands. J. Stat. Phys.; 2018; 173, pp. 1698-1721. [DOI: https://dx.doi.org/10.1007/s10955-018-2151-5]

9. Chen, Z.; de Gier, J.; Hiki, I.; Sasamoto, T. Exact confirmation of 1d nonlinear fluctuating hydrodynamics for a two-species exclusion process. Phys. Rev. Lett.; 2018; 120, 240601. [DOI: https://dx.doi.org/10.1103/PhysRevLett.120.240601]

10. Brevik, I.; Parashar, P.; Shajesh, K.V. Casimir force for magnetodielectric media. Phys. Rev. A; 2018; 98, 032509. [DOI: https://dx.doi.org/10.1103/PhysRevA.98.032509]

11. Xu, F.; Wang, J. Statistical properties of electrochemical capacitance in disordered mesoscopic capacitors. Phys. Rev. B; 2014; 89, 245430. [DOI: https://dx.doi.org/10.1103/PhysRevB.89.245430]

12. Lenzi, E.; Mendes, R.; Rajagopal, A. Green functions based on tsallis nonextensive statistical mechanics: Normalized q-expectation value formulation. Phys. Stat. Mech. Its Appl.; 2000; 286, pp. 503-517. [DOI: https://dx.doi.org/10.1016/S0378-4371(00)00364-2]

13. Lapolla, A.; Vossel, M.; Godec, A. Time- and ensemble-average statistical mechanics of the Gaussian network model. J. Phys. Math. Theor.; 2021; 54, 355601. [DOI: https://dx.doi.org/10.1088/1751-8121/ac16e4]

14. Cornu, F.; Jancovici, B. The electrical double layer: A solvable model. J. Chem. Phys.; 1989; 90, pp. 2444-2452. [DOI: https://dx.doi.org/10.1063/1.455986]

15. Ferrero, A.; Téllez, G. Two-dimensional two-component plasma with adsorbing impurities. J. Stat. Phys.; 2007; 129, pp. 759-786. [DOI: https://dx.doi.org/10.1007/s10955-007-9401-2]

16. Ferrero, A.; Téllez, G. Screening of an electrically charged particle in a two-dimensional two-component plasma at γ=2. J. Stat. Mech. Theory Exp.; 2014; 2014, P11021. [DOI: https://dx.doi.org/10.1088/1742-5468/2014/11/P11021]

17. Joost, J.P.; Schlünzen, N.; Hese, S.; Bonitz, M.; Verdozzi, C.; Schmitteckert, P.; Hopjan, M. Löwdin’s symmetry dilemma within Green functionstheory for the one-dimensional Hubbard model. Contrib. Plasma Phys.; 2021; 62, e202000220.

18. David, G.; Mayboroda, S. Approximation of Green functions and domains with uniformly rectifiable boundaries of all dimensions. Adv. Math.; 2022; 410, 108717. [DOI: https://dx.doi.org/10.1016/j.aim.2022.108717]

19. On the Minimality of Extra Critical Points of Green Functions on Flat Tori. Int. Math. Res. Not.; 2017; 18, pp. 5591-5608.

20. Nomura, Y.; Darmawan, A.S.; Yamaji, Y.; Imada, M. Restricted boltzmann machine learning for solving strongly correlated quantum systems. Phys. Rev. B; 2017; 96, 205152. [DOI: https://dx.doi.org/10.1103/PhysRevB.96.205152]

21. Salazar, D.S.P. Nonequilibrium thermodynamics of restricted boltzmann machines. Phys. Rev. E; 2017; 96, 022131. [DOI: https://dx.doi.org/10.1103/PhysRevE.96.022131]

22. Novoselov, K.S.; Geim, A.K.; Morozov, S.V.; Jiang, D.; Zhang, Y.; Dubonos, S.V.; Grigorieva, I.V.; Firsov, A.A. Electric field effect in atomically thin carbon films. Science; 2004; 306, pp. 666-669. [DOI: https://dx.doi.org/10.1126/science.1102896] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15499015]

23. Fiori, G.; Bonaccorso, F.; Iannaccone, G.; Palacios, T.; Neumaier, D.; Seabaugh, A.; Banerjee, S.; Colombo, L. Electronics based on two-dimensional materials. Nat. Nanotechnol.; 2014; 9, pp. 768-779. [DOI: https://dx.doi.org/10.1038/nnano.2014.207] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25286272]

24. Schwierz, F.; Pezoldt, J.; Granzner, R. Two-dimensional materials and their prospects in transistor electronics. Nanoscale; 2015; 7, pp. 8261-8283. [DOI: https://dx.doi.org/10.1039/C5NR01052G] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25898786]

25. Sterling, R.; Rattanasonti, H.; Weidt, S.; Lake, K.; Srinivasan, P.; Webster, S.; Kraft, M.; Hensinger, W. Fabrication and operation of a two-dimensional ion-trap lattice on a high-voltage microchip. Nat. Commun.; 2014; 5, 3637. [DOI: https://dx.doi.org/10.1038/ncomms4637] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24704758]

26. Flindt, C.; Mortensen, N.A.; Jauho, A.-P. Quantum computing via defect states in two-dimensional antidot lattices. Nano Lett.; 2005; 5, pp. 2515-2518. [DOI: https://dx.doi.org/10.1021/nl0518472]

27. Sugino, F. Super yang-mills theories on the two-dimensional lattice with exact supersymmetry. JHEP; 2004; 403, [DOI: https://dx.doi.org/10.1088/1126-6708/2004/03/067]

28. Saraví, R.G.; Schaposnik, F.; Solomin, J. Path-integral formulation of two-dimensional gauge theories with massless fermions. Nucl. Phys. B; 1981; 185, pp. 239-253. [DOI: https://dx.doi.org/10.1016/0550-3213(81)90375-8]

29. Perera, A.; Urbic, T. Clustering in complex ionic liquids in two dimensions. J. Mol. Liq.; 2018; 265, pp. 307-315. [DOI: https://dx.doi.org/10.1016/j.molliq.2018.05.133]

30. Rañada, M.F.; Santander, M. Superintegrable systems on the two-dimensional sphere s2 and the hyperbolic plane h2. J. Math. Phys.; 1999; 40, pp. 5026-5057. [DOI: https://dx.doi.org/10.1063/1.533014]

31. Kalnins, E.G.; Kress, J.M.; Miller, W. Second-order superintegrable systems in conformally flat spaces. i. two-dimensional classical structure theory. J. Math. Phys.; 2005; 46, 053509. [DOI: https://dx.doi.org/10.1063/1.1897183]

32. Speight, J.M. Static intervortex forces. Phys. Rev. D; 1997; 55, pp. 3830-3835. [DOI: https://dx.doi.org/10.1103/PhysRevD.55.3830]

33. Mitrea, D.; Mitrea, I. On the regularity of green functions in lipschitz domains. Commun. Partial. Differ. Equ.; 2010; 36, pp. 304-327. [DOI: https://dx.doi.org/10.1080/03605302.2010.489629]

34. Truhlar, D.G. Finite difference boundary value method for solving one-dimensional eigenvalue equations. J. Comput. Phys.; 1972; 10, pp. 123-132. [DOI: https://dx.doi.org/10.1016/0021-9991(72)90094-0]

35. Jomaa, Z.; Macaskill, C. The embedded finite difference method for the poisson equation in a domain with an irregular boundary and dirichlet boundary conditions. J. Comput. Phys.; 2005; 202, pp. 488-506. [DOI: https://dx.doi.org/10.1016/j.jcp.2004.07.011]

36. Steger, J.L. Coefficient matrices for implicit finite difference solution of the inviscid fluid conservation law equations. Comput. Methods Appl. Mech. Eng.; 1978; 13, pp. 175-188. [DOI: https://dx.doi.org/10.1016/0045-7825(78)90056-7]

37. Ozis, T.; Aksan, E.; Özdeş, A. A finite element approach for solution of burgers’ equation. Appl. Math. Comput.; 2003; 139, pp. 417-428.

38. Lin, Y.; Xu, C. Finite difference/spectral approximations for the time-fractional diffusion equation. J. Comput. Phys.; 2007; 225, pp. 1533-1552. [DOI: https://dx.doi.org/10.1016/j.jcp.2007.02.001]

39. Izadian, J.; Ranjbar, N.; Jalili, M. The generalized finite difference method for solving elliptic equation on irregular mesh. World Appl. Sci. J.; 2013; 21, pp. 95-100.

40. Jo, G.; Kwak, D. Geometric multigrid algorithms for elliptic interface problems using structured grids. Numer. Algorithms; 2018; 81, pp. 211-235. [DOI: https://dx.doi.org/10.1007/s11075-018-0544-9]

41. Kwak, D.Y.; Kwon, H.J.; Lee, S. Multigrid algorithm for cell centered finite difference on triangular meshes. Appl. Math. Comput.; 1999; 105, pp. 77-85. [DOI: https://dx.doi.org/10.1016/S0096-3003(98)10094-2]

42. Fornberg, B. Generation of finite difference formulas on arbitrarily spaced grids. Math. Comput.; 1988; 51, pp. 699-706. [DOI: https://dx.doi.org/10.1090/S0025-5718-1988-0935077-0]

43. Forsythe, G.; Wasow, W. Finite Difference Methods for Partial Differential Equations: Applied Mathematics Series; Literary Licensing, LLC: Whitefish, MT, USA, 2013.