# Solution to the Modified Helmholtz Equation for Arbitrary Periodic Charge Densities

^{1}Institute for Advanced Simulation, Forschungszentrum Jülich, Jülich, Germany^{2}JARA-CSD, Jülich, Germany^{3}Jülich Supercomputing Centre, Forschungszentrum Jülich, Jülich, Germany^{4}Peter Grünberg Institute, Forschungszentrum Jülich, Jülich, Germany^{5}Physics Department, RWTH-Aachen University, Aachen, Germany

We present a general method for solving the modified Helmholtz equation without shape approximation for an arbitrary periodic charge distribution, whose solution is known as the Yukawa potential or the screened Coulomb potential. The method is an extension of Weinert’s pseudo-charge method [Weinert M, *J Math Phys*, 1981, 22:2433–2439] for solving the Poisson equation for the same class of charge density distributions. The inherent differences between the Poisson and the modified Helmholtz equation are in their respective radial solutions. These are polynomial functions, for the Poisson equation, and modified spherical Bessel functions, for the modified Helmholtz equation. This leads to a definition of a modified pseudo-charge density and modified multipole moments. We have shown that Weinert’s convergence analysis of an absolutely and uniformly convergent Fourier series of the pseudo-charge density is transferred to the modified pseudo-charge density. We conclude by illustrating the algorithmic changes necessary to turn an available implementation of the Poisson solver into a solver for the modified Helmholtz equation.

## 1 Introduction

A variety of problems in condensed matter physics require an efficient solution of the partial differential equation

for a charge density *ρ* in a periodic domain. In real solids, the electronic charge density is a strongly oscillating function in the vicinity of the nuclei of atoms, making a solution in Fourier space, as anticipated by the periodicity, unfeasible due to slow convergence of the Fourier series of the charge density. This equation is frequently referred to as the modified Helmholtz equation or the Yukawa equation. The latter name derives from the Yukawa potential [1], *e.g.* in physics, chemistry, and biology, the Yukawa potential is also known as the screened Coulomb potential. It typically emerges in cases when a many-body system of charged particles is treated in terms of an effective single-particle theory applying a mean-field approximation. Then the many particles contribute to an effective screening of a Coulomb interaction generated by the single, representative charged particle when treated in linear response theory.

The relation between the bare Coulomb potential on the one hand and the screened Coulomb potential or the induced screening charge on the other hand is referred to as the dielectric constant or susceptibility, respectively. Depending on the context, such relations appear in the Debye–Hückel theory [2] in the form of a linearization of the Poisson–Boltzmann equation, where the Poisson equation describing the electrostatics of charged particles is a function of the charge density distribution obeying a Boltzmann statistics. Another example is the Thomas–Fermi model [3, 4] of the dielectric constant in metals, which describes the screening potential due to the linearized change of the electron distribution described by the Fermi–Dirac distribution with respect to the spatial variations of the electrostatic potential. In these theories, the constant *λ* represents the inverse of a typical length scale over which an individual charged particle exerts a notable effect.

The Thomas-Fermi theory can be regarded as a precursor of the density functional theory [5] (DFT). The latter is the most important theory and methodology for the modeling and simulation of material properties of a crystalline solid based on the quantum mechanical treatment of many electron systems. In addition, the Thomas-Fermi theory provides a rough but fast approximation of the common density functionals, which relate the electron density to the effective Kohn–Sham potentials [6]. Such a scheme makes the solution of Eq. 1 particularly valuable. For instance it could be used to obtain a good starting potential for the iterative solution of the Schrödinger-like Kohn–Sham equations, where the nuclear charge is included in *ρ*. Other examples are the attainment of an efficient approximate solution of the dielectric function, or the implementation of a hybrid functional [7] to DFT using the Yukawa screening of the Hartree–Fock exchange. In both cases the charge density *ρ* in Eq. 1 is replaced by an overlap charge density [8] obtained as a product of wave functions associated with different quantum numbers. Since the Thomas–Fermi model approximates the description of the response behavior of delocalized electrons in solids very well, Eq. 1 is frequently used to determine a preconditioner for the acceleration of the self-consistent solution of the DFT equation. So far the Kerker preconditioner [9] is applied to density functional methods with Fourier transformable charge density. By solving Eq. 1 for general periodic densities, the Kerker preconditioner [9] can be extended to density functional methods of general densities [10, 11].

Although most electronic structure methods implementing DFT applied to solid-state materials systems make explicit use of the underlying periodicity of the crystalline lattice, a straightforward solution of Eq. 1 using Fourier transformation techniques is in general not possible due to the strongly oscillating charge density close to the nuclei. This problem is well discussed for the solution of the Poisson equation,

In a seminal work, Weinert [12] proposed an elegant and numerically efficient solution of the Poisson equation for periodic charges and corresponding electrostatic potentials without shape approximation. Weinert’s solution, to which we refer here as Weinert’s pseudo-charge method, is implemented (in several variants) in most full-potential all-electron DFT methods, such as the augmented spherical wave (ASW) method [13], the Korringa–Kohn–Rostoker Green function (KKR-GF) method [14], and the full-potential linearized augmented planewave (FLAPW) method [15], just to name a few.

Typical to these all-electron DFT-methods is the domain decomposition into atomic spheres around the atoms and an interstitial region in-between. Weinert’s pseudo-charge method is based on the observation that the relation between the charge density inside a sphere and its multipole expansion outside the sphere is not unique. A smooth Fourier transformable pseudo-charge density with the same multipole moments as the true density is constructed. The latter provides the true potential through Fourier transformation of the Poisson equation and a subsequent solution of a Dirichlet boundary value problem on the sphere boundary.

In this article, we extend Weinert’s pseudo-charge method to the modified Helmholtz Eq. 1 for values of

As a matter of choice, and motivated by the original work of Weinert [12], we demonstrate this extension explicitly for the FLAPW method [15] as implemented in the Fleur code [16]. We provide i) complete derivations of the modified multipole expansion of the charge density inside the atomic spheres using a Green function method, as well as the interstitial charge density’s modified multipole moments in the atomic spheres using some Bessel function integration properties, which yields the coefficients of the pseudo-charge density, and ii) a mathematical analysis of the convergence of the pseudo-charge density’s Fourier series. Our results confirm all expressions for all the quantities provided by Tran et al. [7] relevant in an actual computation. Beyond Tran et al. [7] we provide detailed derivations and justify the use of identities such as Proposition 1 for the modified spherical Bessel function

This paper is organized as follows: In Section 2.1, we introduce the muffin-tin and interstitial region typical of the FLAPW method and the corresponding domain decomposition for the charge density and the potential. We summarize the main statements of Weinert’s pseudo-charge method and give a definition of the pseudo-charge density. Since we know the true charge density inside the muffin-tin spheres and with the assumption that we would know the interstitial potential, we construct in Section 2.2 the Yukawa potential inside the sphere by solving the Dirichlet boundary value problem. We develop two radial Green functions that are products of two linearly independent fundamental set solutions of the homogeneous radial modified Helmholtz equation. These Green functions are set apart by the boundary conditions they fulfill either at the muffin-tin sphere or in free-space. In Section 2.2.1, the radial free-space Green function is used to define the modified multipole expansion of the Yukawa potential. In Section 2.3, we construct a pseudo-charge density in reciprocal space consistent with the modified Helmholtz equation by making use of the definition of the modified multiple moments put forward in Section 2.2.1. We obtain the Yukawa potential for the interstitial region by solving the modified Helmholtz equation in Fourier space for the pseudo-charge density—the solution is a simple algebraic expression. This is followed by an analysis of the convergence properties of the Fourier series of the pseudo-charge density. The entire algorithm that solves the modified Helmholtz equation is summarized in Section 2.4 together with the minimal modifications necessary to change Weinert’s original algorithm. The conclusions and the outlook are presented in Section 3.

## 2 Yukawa Potential for a Muffin-Tin Decomposition of a 3D-Periodic Domain

### 2.1 Weinert’s Pseudo-Charge Method

In order to deal with the *α*—the union of those is called the muffin-tin (MT) region—and the interstitial region (I) between the atoms. In FLAPW both charge densities

and potentials

are represented in plane waves *Y*_{L}, of degree *m*, where *L* is defined as ** K**, with length

Weinert’s pseudo-charge method for the Poisson equation is based on the crucial observation that several charge densities *ρ* inside a sphere

and thus, the same potential

outside the sphere. Here,

• It has the same multipole moments

• It is equal to the true charge density

• It has a fast convergent Fourier expansion.

The Fourier Components of

while

where *G* is a Green function associated with the solution of the Poisson equation, *R*_{α}, for simplicity we drop the index α in the Green function and in related quantities. The muffin-tin potential

where *i.e.* including the muffin-tin spheres, such that

can also be evaluated for *ρ* can also be written as

where

are charge densities localized in the atomic spheres

The approach above is generally the same for the Yukawa potential, *i.e.,* for

### 2.2 Muffin-Tin Yukawa Potential

Assume the interstitial potential is obtained, then the Green function method is used to determine the screened Coulomb potential

The solution is divided into three steps, which we describe in the following. The derivation focuses on the construction of the Green function and its application. For simplicity, we leave out the index *λ* in Green functions and potentials in this subsection. The solution is given in terms of radial functions

*Step 1*. We solve the homogeneous modified Helmholtz equation,

and its fundamental set of solutions are for each *r* and the second is the irregular solution that goes to infinity for

since

*Step 2*. A function

subject to the Dirichlet boundary condition

where

The radial Green function takes the form of the product of the two linearly independent solutions with the proper boundary conditions,

where

Since the Wronskian *W* is linear, the addition of *c* = const) onto

Therefore,

*Step 3*. Considering the standard expression of Dirac’s delta function separated according to the radial and angular coordinates

the three-dimensional (3D) Green function

is expanded in the form

and the solution to the inhomogeneous Eq. 12 is given by Eq. 7. For the derivation of the 3D Green function and the 3D inhomogeneous solution

In order to obtain the boundary contribution to the muffin-tin potential, we evaluate the boundary integral in Eq. 7 by expanding the interstitial potential

using the plane-wave expansion

where

Furthermore, the normal derivative of *G* on the sphere boundary is

Since

We recall that

yielding

With this and the knowledge of the interstitial potential at the sphere boundary from Eq. 31, the boundary contribution to the muffin-tin potential becomes

and the radial part of the spherical harmonics expansion of the total potential,

Due to the kink of

##### 2.2.1 Modified Multipole Expansion

In the same way we obtain the radial representation of the free-space Green function, well-known as the Yukawa potential for a Dirac test charge at

The modified spherical Bessel function *i.e*.,

which can be expressed analogously to the Coulomb potential in terms of the modified multipole expansion

with the modified multipole moments

An analogous definition holds true for the modified multiple moments

The summation of

### 2.3 Interstitial Yukawa Potential

Suppose we had found a Fourier transformable pseudo-charge density

In the previous Section 2.2 the interstitial Yukawa potential is used as boundary values for the Yukawa potential in the atomic spheres.

This subsection is concerned with the construction of the Fourier transformable pseudo-charge density *a priori* unknown pseudo-charge density,

##### 2.3.1 Modified Multipole Moments of Interstitial Charge Density Extended Into Sphere

The determination of the modified multiple moments of the interstitial charge density is in principle the same as in Ref. [12], but since the modified multipole moments are different from the known multipole moments for the Coulomb potential, we go through this step of deriving

The modified multipole moments of

For

If *δ*) and so the integral becomes

Both equations above can be derived by partial integration in two different ways and applying the identities in Ref. [18],

for

for both

##### 2.3.2 Construction of Pseudo-Charge Density

We construct the pseudo-charge density by following the ansatz of Weinert,

in which the radial dependence of the charge density is expressed in terms of a polynomial expansion up to degree

and *n* is yet to be determined. As will become apparent later when we derive the coefficients

it follows from the ansatz Eq. 51

The Fourier transform of this expression is then given by

where

In the same way we derive the coefficients

and thus,

Since

and setting

For

which is only different from 0 for

Due to the condition that the pseudo-charge density has the correct modified multipole moments, the modified multipole moments used for the actual computation of the Fourier coefficients are the ones of the true localized charge density

##### 2.3.3 Smoothness of the Pseudo-Charge Density and Convergence of Its Fourier Series

In Section 2.3.2, we have set *n* yet. Here we motivate our choices and finally determine a proper *n*.

With our choices for

itself and all its first

Note that this includes localization of

The smoothness of the pseudo-charge density is connected to the convergence properties of its Fourier series. Applying the Riemann–Lebesgue lemma for Fourier series to a *K* and the convergence becomes the better the larger *n* is. If we choose *n* too large, however, we are left with the small-*K* Fourier coefficients only, and thus the Fourier series is unbalanced, in the sense that smaller coefficients have a larger weight. So, ideally, our choice of *n* is guided by the cut-off of the Fourier series.

For an explicit rule on how to choose *n*, Weinert [12] discussed the factors

His reasoning is based on the zeros of the *K*-dependent function *K*, however, our factor *n* and angular degree *n*, and (ii) for fixed *n*, the largest contribution to the Fourier series comes from ** K** vectors, we would like to reduce the oscillations mentioned in (i) by choosing a large

*n*. On the other hand, due to the cut-off of the Fourier series at some

• Choose

• Then

**FIGURE 1**. The factor

Note that in this method *n* in such a way that *ν* is *de facto* not depending on ** K** vectors and the muffin-tin radii usually do not change over the course of the self-consistent-field iteration, the terms in Eq. 63 depending on

*ν*need to be computed just once.

### 2.4 Algorithm: Construction of Yukawa Potential

Algorithm 1 summarizes the construction of the Yukawa potential derived in this paper.

Algorithm 1Bulk-Case Yukawa Potential

**Input:** charge density *ρ*, integer *ν* chosen as described on page 8 and preconditioning parameter *λ*.

**Output:** Yukawa potential

**Pseudo-Charge Density**

1: Modified multipole moments

2: Modified multipole moments

3: Modified Multipole Moments

4: Sphere-localized part

5: Pseudo-charge density **Interstitial potential**

6: Interstitial potential **Muffin-tin potential**

7: Boundary terms

8: Radial parts

In the case that Weinert’s method is available as an implemented algorithm, then only relatively few changes are necessary to make it available for the solution of the modified Helmholtz equation. The changes to be made in practice are limited to the following: The slightly different radial behavior of the Green function leads to small changes in the multipole moments of the interstitial and muffin-tin charge densities in each sphere, *ν* in the formula for the pseudo-charge density’s Fourier components, which determines the convergence of the Fourier series, is chosen exactly the same as in Weinert’s original method. The interstitial potential Eq. 44 undergoes changes indirectly through the pseudo-charge density and directly by the prefactor ** K** =

**0**-term is well-defined, it is not set to a constant as in the original method. The muffin-tin potential is affected only in its radial dependence in both the boundary and the source contribution of Eq. 2.2. Basically, the polynomials

## 3 Conclusion and Outlook

We have presented a general method for solving the modified Helmholtz equation for a 3D-periodic system of charge densities not restricted by any shape approximation of three-dimensional volume. The three-dimensional domain is typically decomposed into non-overlapping atom-centered spheres (the muffin-tin region) and the space between these spheres (the interstitial region). The solution, the Yukawa potential, suffices 3D-periodic boundary conditions. Since the Yukawa differential equation is similar to the Poisson equation, we leveraged our derivations on the work of Weinert [12]. Our work can be considered an extension of Weinert’s work determining instead of the bare Coulomb potential with zero screening, the screened Coulomb potential with finite screening length

Considering that Weinert’s pseudo-charge method has become the standard method for calculating the electrostatic potential without shape approximation in all-electron band structure methods for applications of periodic solids, and since we have extended it to a modified pseudo-charge method with only minor modifications involving some radial integrals, this now allows to treat the screened Coulomb potentials without shape approximation described by the modified Helmholtz equation with all-electron methods. The screened Coulomb or Yukawa potential typically occurs in single particle or mean-field theories to the problem of many charged particles, where all charged particles contribute effectively to the screening of the bare Coulomb potential.

## Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

## Author Contributions

MW did all the derivations, provided the proof and the figure. EN and DW supervised the work. SB conceptualized the work. All authors contributed through discussions and to the writing of the paper.

## Funding

This work has been supported by a JARA-HPC seed-fund project and by the MaX Center of Excellence funded by the EU through the H2020-INFRAEDI-2018-1 (Grant No. 824143). The authors gratefully acknowledge the computing time granted by the JARA-HPC Vergabegremium and provided on the JARA-HPC partition’s part of the supercomputer JURECA 20 at Forschungszentrum Jülich, as well as the funds received for open access publication fees from the Forschungszentrum Jülich GmbH.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

We would like to thank Gregor Michalicek, Jan Winkelmann, and Rudolf Zeller for fruitful discussions and help with conceptual and computational matters. This manuscript has been released as a pre-print at Hinzen et al. [21].

## Appendix

Proposition 1. *Let**,**and**be either the spherical Bessel function of the first kind,**, or the modified spherical Bessel function of the first kind,**. Then*

*holds for all**and all*

Proof by mathematical induction on *n*. The base case

The induction step uses the property for *n* and

be the left-hand side of Eq. 1. Then

where we used the induction hypothesis in the last equation, and thus the proposition follows. ∎Proposition 2. For nonnegative λ and R, and

*holds, with equality if and only if**or*

where

and it remains to show that

The first term of the sum for *λ* and *R*. The other summed terms for *λ* and *R* are positive, and zero, if one of the two variables is zero. ∎

## References

1. Yukawa H. On the interaction of elementary particles. I. *Proc Physico-Math Soc Jpn* (1935) 33:48–57. doi:10.11429/ppmsj1919.17.0_48

2. Hückel E, Debye P. The theory of electrolytes. I. lowering of freezing point and related phenomena. *Phys Z* (1923) 24:185–206.

3. Thomas L. The calculation of atomic fields. *Math Proc Cambridge Philos Soc* (1927) 23:542–8. doi:10.1017/S0305004100011683

4. Fermi E. Eine statistische Methode zur Bestimmung einiger Eigenschaften des Atoms und ihre Anwendung auf die Theorie des periodischen Systems der Elemente. *Zeitschrift für Physik* (1927) 48:73–9. doi:10.1007/BF01351576

5. Hohenberg P, Kohn W. Inhomogeneous electron gas. *Phys Rev* (1964) 136:B864–71. doi:10.1103/PhysRev.136.B864

7. Tran F, Blaha P. Implementation of screened hybrid functionals based on the Yukawa potential within the LAPW basis set. *Phys Rev B* (2011) 83:235118. doi:10.1103/PhysRevB.83.235118

8. Massidda S, Posternak M, Baldereschi A. Hartree-Fock lapw approach to the electronic properties of periodic systems. *Phys Rev B* (1993) 48:5058–68. doi:10.1103/PhysRevB.48.5058

9. Kerker G. Efficient iteration scheme for self-consistent pseudopotential calculations. *Phys Rev B* (1981) 23:3082.

10. Winkelmann M, Di Napoli E, Wortmann D, Blügel S. Kerker mixing scheme for self-consistent muffin-tin based all-electron electronic structure calculations. *Phys Rev B* (2020) 102:195138. doi:10.1103/PhysRevB.102.195138

11. Kim J, Gulans A, Draxl C. Robust mixing in self-consistent linearized augmented planewave calculations. *Electron Struct* (2020) 2:037001. doi:10.1088/2516-1075/ababde

12. Weinert M. Solution of Poisson’s equation: beyond Ewald-type methods. *J Math Phys* (1981) 22:2433–9. doi:10.1063/1.524800

13. Eyert V. *The plane-wave based full-potential ASW method*. Berlin, Heidelberg, Germany: Springer (2013). p. 113–72. doi:10.1007/978-3-642-25864-0˙4

14. Drittler B. *KKR-Greensche Funktionsmethode für das volle Zellpotential*. [PhD thesis]. Aachen (Germany): RWTH Aachen (1991).

15. Wimmer E, Krakauer H, Weinert M, Freeman AJ. Full-potential self-consistent linearized-augmented-plane-wave method for calculating the electronic structure of molecules and surfaces: O_{2} molecule. *Phys Rev B* (1981) 24:864–75. doi:10.1103/PhysRevB.24.864

16.FLEUR. *The Jülich FLAPW code family*. Available from: http://www.flapw.de/ (2020).

17. Abramowitz M, Stegun IA. “*Handbook of mathematical functions with formulas, graphs, and mathematical tables,” in Ninth dover printing, tenth gpo printing*. New York, NY: Dover (1964).

18. Arfken G, Weber H, Harris F. *Mathematical methods for physicists: a comprehensive guide*. Amsterdam, Netherlands: Elsevier Science (2013).

Keywords: partial differential equations, density functional theory, electronic structure methods, Green functions technique, materials science, electrostatics, Fourier analysis, muffin-tin approximation

Citation: Winkelmann M, Di Napoli E, Wortmann D and Blügel S (2021) Solution to the Modified Helmholtz Equation for Arbitrary Periodic Charge Densities. *Front. Phys.* **8**:618142. doi: 10.3389/fphy.2020.618142

Received: 16 October 2020; Accepted: 16 December 2020;

Published: 11 March 2021.

Edited by:

Manuel Asorey, University of Zaragoza, SpainReviewed by:

Kazuharu Bamba, Fukushima University, JapanBabak Shiri, Neijiang Normal University, China

Copyright © 2021 Winkelmann, Di Napoli, Wortmann and Blügel. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Miriam Winkelmann, m.winkelmann@fz-juelich.de