Front. Phys. Frontiers in Physics Front. Phys. 2296-424X Frontiers Media S.A. 618142 10.3389/fphy.2020.618142 Physics Original Research Solution to the Modified Helmholtz Equation for Arbitrary Periodic Charge Densities Winkelmann et al. Solution to Modified Helmholtz Equation Winkelmann Miriam 1 2 3 4 5 * Di Napoli Edoardo 1 2 3 Wortmann Daniel 1 2 4 Blügel Stefan 1 2 4 Institute for Advanced Simulation, Forschungszentrum Jülich, Jülich, Germany JARA-CSD, Jülich, Germany Jülich Supercomputing Centre, Forschungszentrum Jülich, Jülich, Germany Peter Grünberg Institute, Forschungszentrum Jülich, Jülich, Germany Physics Department, RWTH-Aachen University, Aachen, Germany

Edited by: Manuel Asorey, University of Zaragoza, Spain

Reviewed by: Kazuharu Bamba, Fukushima University, Japan

Babak Shiri, Neijiang Normal University, China

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

This article was submitted to Mathematical and Statistical Physics, a section of the journal Frontiers in Physics

 11 03 2021 2020 8 618142 16 10 2020 16 12 2020 Copyright © 2021 Winkelmann, Di Napoli, Wortmann and Blügel. 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.

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.

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

A variety of problems in condensed matter physics require an efficient solution of the partial differential equation ( Δ − λ 2 ) V λ = − 4 π ρ   , 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], V λ ∝ exp ( − λ r ) / r , in nuclear physics, which is the underlying free-space Green function of Eq. 1. In the field of condensed matter, 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, Δ V = − 4 π ρ , a limit of the modified Helmholtz equation for λ = 0 .

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 λ > 0 in order to determine the Yukawa potential V λ for general periodic charge densities without shape approximation. We formulate the new method for general charge densities, including continuous charge densities as for electron densities, discrete charge densities as for nuclear charges or more abstract densities that arise of products of wave functions. Such an approach is consistent with the real-space representation of the charge density and potential in all-electron methods.

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 i ℓ and the transferability of Weinert’s convergence analysis to the case λ > 0 . We also discuss and provide further insight into the convergence properties of the Fourier series, and in addition restructure and simplify Weinert’s approach, hopefully making it more comprehensible. We point out the algorithmic changes required to extend the solution of the Poisson equation to a solution of Eq. 1, which can then be straightforwardly transferred to other all-electron full-potential band structure methods.

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 1 / r singularities of the Coulomb potential due to the point-like charge of the nucleus and the associated rapid oscillations of the charge density in the vicinity of the singularity, in all-electron electronic structure methods the space is typically partitioned into muffin-tin spheres B R α ( τ α ) of radius R α centered around the atoms α—the union of those is called the muffin-tin (MT) region—and the interstitial region (I) between the atoms. In FLAPW both charge densities ρ ( r ) = { ∑ K ρ I ( K )   e i K ⋅ r     r ∈ I ∑ L ρ L α ( r α )   Y L ( r ^ α ) r = τ α + r α ∈ B R α ( τ α ) and potentials V ( r ) = { ∑ K V I ( K )   e i K ⋅ r   r ∈ I ∑ L V L α ( r α )   Y L ( r ^ α )   r = τ α + r α ∈ B R α ( τ α ) are represented in plane waves e i K ⋅ r , where K defines the reciprocal lattice vector dual to the lattice vectors defining the periodic domain, and in spherical harmonics, Y L , of degree ℓ and order m, where L is defined as L : = ( ℓ , m ) . r α ≤ R α is the length of the vector r α = r − τ α , measured from the center of the atom α placed at position τ α in the periodic domain, with r ^ α = r − τ α / | r − τ α | its unit vector. The precision of the representation is determined by the cut-off parameters K max for the wave vectors, K , with length K ≤ K max , and ℓ max for the degree ℓ in the angular-momentum expansion. ℓ max sets also a natural cut-off of the angular-momentum expansions of all other charge densities or multipole moments throughout the paper.

Weinert’s pseudo-charge method for the Poisson equation is based on the crucial observation that several charge densities ρ inside a sphere B R α ( τ α ) can generate the same multipole moments q L α [ ρ ] = ∫ B R α ( 0 ) ρ ( r α + τ α ) r α ℓ Y L * ( r ^ α ) d r α and thus, the same potential V I [ ρ ] ( r ) = ∑ L 4 π 2 ℓ + 1   q L α [ ρ ] 1 r α ℓ + 1 Y L ( r ^ α ) outside the sphere. Here, Y L * denotes the complex conjugate of Y L . The pseudo-charge density, ρ ˜ , defined by Weinert in Ref. [12] is such a charge density. It fulfills the following three conditions:

• It has the same multipole moments q L α [ ρ ˜ ] = q L α [ ρ ] in every sphere B R α ( τ α ) .

• It is equal to the true charge density ρ I in the interstitial region.

• It has a fast convergent Fourier expansion.

The Fourier Components of V I are then simply V I ( K ) = 4 π K 2 ρ ˜ ( K )   for   K ≠ 0 , while V I ( 0 ) will be set to a constant. Once the interstitial potential V I has been calculated, the muffin-tin potential can be obtained by solving the Dirichlet boundary value problem on the sphere V α ( r α + τ α ) = V S α ( r α + τ α ) + V B α ( r α + τ α ) = ∫ B R α ( 0 ) G ( r α , r α ′ ) ρ ( r α ′ + τ α ) d r α ′ − R α 2 4 π ∫ ∂ B R α ( 0 ) V I ( r α ′ + τ α ) ∂ G ∂ n ′ ( r α , r α ′ ) d ω ′ , where G is a Green function associated with the solution of the Poisson equation, d ω = sin θ d θ d ϕ denotes the solid angle element and r = r α + τ α ∈ B R α ( τ α ) . Although the Green function depends on the muffin-tin radius R α , for simplicity we drop the index α in the Green function and in related quantities. The muffin-tin potential V α is fed by two terms, a source term V S α due to the charge density distribution inside the sphere and a boundary term V B α due to the interstitial potential at the boundary of the sphere. The Fourier coefficients of the pseudo-charge density basically have the form ρ ˜ ( K ) = ρ I ( K ) + ∑ α ρ ¯ α ( K )   , where ρ ¯ α is a Fourier transformable pseudo-charge density inside the muffin-tin sphere. The idea behind this is the following: if the domain of definition of ρ I is formally expanded to the full space, i.e. including the muffin-tin spheres, such that ρ I ( r ) = ∑ K ρ I ( K ) e i K ⋅ r can also be evaluated for r ∈ ∪ α B R α ( τ α ) , then the true charge ρ can also be written as ρ = ρ I + ∑ α ρ ˘ α   , where ρ ˘ α : = { 0           in   I ρ α − ρ I   in   B R α ( τ α ) are charge densities localized in the atomic spheres B R α ( τ α ) . If these localized densities are now substituted by other localized densities ρ ¯ α , then the charge density is still correct in I.

The approach above is generally the same for the Yukawa potential, i.e., for λ > 0 , only Eqs. 4–6 and the pseudo-charge density’s Fourier coefficients are slightly different. We now continue with the derivation of the muffin-tin and interstitial potentials in the case λ > 0 .

 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 V λ α inside the muffin-tin sphere B R α ( τ α ) with centre τ α through the solution of the boundary value problem ( Δ − λ 2 ) V λ α = − 4 π ρ in  B R α ( τ α ) V λ α = V λ I on ∂ B R α ( τ α ) . 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 V L α ( r α ) , the expansion coefficients to the spherical harmonics expansion inside the sphere (see Eq. 3).

Step 1. We solve the homogeneous modified Helmholtz equation, ( Δ − λ 2 ) U = 0 , in spherical coordinates. Following the solution [17] of the Laplace equation, Δ ψ = 0 , in spherical coordinates ( r , r ^ ) , the homogeneous potential can be factorized into products of radial functions u ℓ ( r ) and angular functions Y L ( r ^ ) , U ( r , r ^ ) = ∑ L u ℓ ( r ) Y L ( r ^ ) . The term − λ 2 U in the homogeneous modified Helmholtz equation only has an effect on the radial solution. The spherical harmonics, Y L , are the eigensolutions of the angular part of the Laplace equation with eigenvalues ℓ ( ℓ + 1 ) . The radial part of the homogeneous modified Helmholtz equation is known as the modified spherical Bessel differential equation [17, 18], d 2 u ℓ ( r ) d r 2 + 2 r d u ℓ ( r ) d r − ( ℓ ( ℓ + 1 ) r 2 + λ 2 ) u ℓ ( r ) = 0 , and its fundamental set of solutions are for each ℓ the two modified spherical Bessel functions [17, 18] i ℓ ( λ r ) and k ℓ ( λ r ) , the first of which is the regular solution well-defined at the origin, but grows fast with growing radius r and the second is the irregular solution that goes to infinity for r → 0 . To realize the proper boundary condition for the radial Green functions two conditions have to be fulfilled: The first solution, u ℓ 1 , must be finite at r = 0 . We conclude that u ℓ 1 ( r ) = i ℓ ( λ r )   , since k ℓ ( λ r ) → ∞ for r → 0 . The second solution, u ℓ 2 , must be 0 at r = R α . This is achieved by a linear combination of the two modified spherical Bessel functions, u ℓ 2 ( r ) = k ℓ ( λ r ) − i ℓ ( λ r ) k ℓ ( λ R α ) i ℓ ( λ R α )   .

Step 2. A function G ℓ ∈ C 0 ( [ 0 , R α ] × [ 0 , R α ] ) ∩ ​ C 2 ( [ 0 , R α ] × [ 0 , R α ] \ { ( r , r ) | r ∈ [ 0 , R α ] } ) is called a radial Green function, if it is the solution to ( Δ r − λ 2 ) G ℓ ( r , r ′ ) = − 4 π r 2 δ ( r − r ′ ) subject to the Dirichlet boundary condition G ℓ ( r , r ′ ) = 0   ,     if   r = R α   or   r ′ = R α   , where ( Δ r − λ 2 ) is the linear radial differential operator in Eq. 14 and δ denotes the radial Dirac delta function, for which ∫ ​ a b δ ( r − r ′ ) f ( r ′ ) d r ′ = { f ( r ) , if  r ∈ [ a , b ] 0 , otherwise .   The radial Green function takes the form of the product of the two linearly independent solutions with the proper boundary conditions, G ℓ ( r , r ′ ) = C u ℓ 1 ( r < ) u ℓ 2 ( r > )   , where r < = min ( r , r ′ ) , r > = max ( r , r ′ ) and C = − 4 π r ′ 2 W − 1 [ u ℓ 1 ( r ′ ) , u ℓ 2 ( r ′ ) ] . Since the Wronskian W is linear, the addition of c u ℓ 1 (c = const) onto k ℓ ( λ r ) to suffice the boundary condition u ℓ 2 ( R α ) = 0 has no influence on the Wronskian: W ( u ℓ 1 , c u ℓ 1 ) = 0 , with W ( i ℓ , k ℓ ) remaining. To calculate the Wronskian of i ℓ ( λ r ′ ) and k ℓ ( λ r ′ ) one can either Taylor expand the two functions or simply take the limiting values for r → 0 or the asymptotic values for r → ∞ to find W [ u ℓ 1 ( r ′ ) , u ℓ 2 ( r ′ ) ] = u ℓ 1 ( r ′ ) d u ℓ 2 ( r ′ ) d r ′ − d u ℓ 1 ( r ′ ) d r ′ u ℓ 2 ( r ′ ) = − 1 λ r ′ 2   . Therefore, C = 4 π λ , and the radial Green function finally reads G ℓ ( r , r ′ ) = 4 π λ   i ℓ ( λ r < ) k ℓ ( λ r > ) − 4 π λ   i ℓ ( λ r ) i ℓ ( λ r ′ ) k ℓ ( λ R α ) i ℓ ( λ R α )   .

Step 3. Considering the standard expression of Dirac’s delta function separated according to the radial and angular coordinates δ ( r − r ′ ) = 1 r 2 δ ( r − r ′ ) δ ( r ^ − r ^ ′ ) = 1 r 2 δ ( r − r ′ ) ∑ L Y L * ( r ^ ′ ) Y L ( r ^ )   , the three-dimensional (3D) Green function G ( r , r ′ ) solving ( Δ − λ 2 ) G ( r , r ′ ) = − 4 π δ ( r − r ′ ) in  B R α ( 0 ) , G ( r , r ′ ) = 0 on  ∂ B R α ( 0 )   , is expanded in the form G ( r , r ′ ) = ∑ L G ℓ ( r , r ′ ) Y L * ( r ^ ′ ) Y L ( r ^ ) 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 V α = V S α + V B α Eq. 7 we refer the reader to Ref. [19]. Both the integral over the sphere B R α ( 0 ) and the boundary integral simplify by exploiting the orthonormality relation of the spherical harmonics, ∫ ∂ B 1 ( 0 ) Y L * ( r ^ ) Y L ′ ( r ^ ) d ω = δ L L ′ . The integral over the sphere B R α ( 0 ) provides the source contribution to the muffin-tin potential V S α ( r α + τ α ) = ∫ B R α ( 0 ) G ( r α , r α ′ ) ρ ( r α ′ + τ α ) d r α ′ = ∑ L [ ∫ ​ 0 R α G ℓ ( r α , r α ′ ) ρ L α ( r α ′ ) r α ′ 2 d r α ′ ] Y L ( r ^ α ) . 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 V I ( K ) (see Section 2.3) on the sphere boundaries ∂ B R α ( τ α ) ∋ r ′ in spherical coordinates V I ( r α ′ + τ α ) = ∑ K V I ( K ) e i K ⋅ τ α e i K ⋅ r α ′ = ∑ L V L I ( R α ; τ α ) Y L ( r ^ α ′ ) using the plane-wave expansion e i K ⋅ r = ∑ L 4 π i ℓ j ℓ ( K r ) Y L * ( K ^ ) Y L ( r ^ )   , where V L I ( R α ; τ α ) = 4 π i ℓ ∑ K V I ( K ) e i K ⋅ τ α j ℓ ( K R α ) Y L * ( K ^ )   . Furthermore, the normal derivative of G on the sphere boundary is ∂ G ∂ n ′ ( r α , r α ′ ) = ∂ G ( r α , r α ′ ) ∂ r α ′ | r α ′ = R α = ∑ L ∂ G ℓ ( r α , r α ′ ) ∂ r α ′ | r α ′ = R α Y L * ( r ^ α ′ ) Y L ( r ^ α )   . Since r α < R α = r α ′ and since G ℓ takes the form Eq. 20, we obtain ∂ G ℓ ( r α , r α ′ ) ∂ r α ′ | r α ′ = R α = 4 π λ   u ℓ 1 ( r α ) u ℓ 2 ′ ( R α )   . We recall that u ℓ 2 ( R α ) = 0 and reuse Eq. 22 to obtain u ℓ 2 ′ ( R α ) = − 1 λ R α 2 i ℓ ( λ R α )   , yielding ∂ G ℓ ( r α , r α ′ ) ∂ r α ′ | r α ′ = R α = − 4 π R α 2 i ℓ ( λ r α ) i ℓ ( λ R α )   . 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 V B α ( r α + τ α ) = − R α 2 4 π ∫ ∂ B R α ( 0 ) V I ( r α ′ + τ α ) ∂ G ∂ n ′ ( r α , r α ′ ) d ω ′ = ∑ L V L I ( R α ; τ α ) i ℓ ( λ r α ) i ℓ ( λ R α ) Y L ( r ^ α ) and the radial part of the spherical harmonics expansion of the total potential, V S α + V B α , in the sphere B R α ( 0 ) becomes V L α ( r α ) = ∫ ​ 0 R α G ℓ α ( r α , r α ′ ) ρ L α ( r α ′ ) r α ′ 2 d r α ′ + V L I ( R α ; τ α ) i ℓ ( λ r α ) i ℓ ( λ R α )   . Due to the kink of G ℓ at r α = r α ′ , for practical calculations the integral is split in a part where r α ′ < r α , a part where r α ′ > r α and a third part where the integrand is symmetric in r α and r α ′ : V L α ( r α ) = 4 π λ ( [ ∫ ​ 0 r α ρ L α ( r α ′ ) i ℓ ( λ r α ′ ) r α ′ 2 d r α ′ ] k ℓ ( λ r α ) + [ ∫ ​ r α R α ρ L α ( r α ′ ) k ℓ ( λ r α ′ ) r α ′ 2 d r α ′ ] i ℓ ( λ r α ) − [ ∫ ​ 0 R α ρ L α ( r α ′ ) i ℓ ( λ r α ′ ) r α ′ 2 d r α ′ ] i ℓ ( λ r α ) k ℓ ( λ R α ) i ℓ ( λ R α ) ) + V L I ( R α ; τ α ) i ℓ ( λ r α ) i ℓ ( λ R α )   .

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 r ′ , e − λ | r − r ′ | | r − r ′ | = 4 π λ ∑ L i ℓ ( λ r < ) k ℓ ( λ r > ) Y L * ( r ^ ′ ) Y L ( r ^ )   . The modified spherical Bessel function k ℓ already contains the proper boundary condition for r → ∞ . A charge density ρ ˘ α localized in a sphere B R α ( τ α ) embedded in free space, produces a Yukawa potential outside the sphere, i.e., r = r α + τ α ∉ B R α ( τ α ) , V I [ ρ ˘ α ] ( r ) = ∫ ℝ 3 G ( r α , r α ′ )   ρ ˘ α ( r α ′ + τ α ) d r α ′   , which can be expressed analogously to the Coulomb potential in terms of the modified multipole expansion V I [ ρ ˘ α ] ( r ) = ∑ L 4 π λ ℓ + 1 ( 2 ℓ + 1 ) ! ! q L α [ ρ ˘ α ]   k ℓ ( λ r α )   Y L ( r ^ α )   , with the modified multipole moments q L α [ ρ ˘ α ] = ( 2 ℓ + 1 ) ! ! λ ℓ ∫ B R α ( 0 ) ρ ˘ α ( r α + τ α ) i ℓ ( λ r α ) Y L * ( r ^ α ) d r α   . An analogous definition holds true for the modified multiple moments q L α [ ρ α ] of the true charge ρ α in the sphere. With the standard expansion of the charge density inside the sphere into spherical harmonics Eq. 3, ρ α ( r α + τ α ) = ∑ L ρ L α ( r α )   Y L ( r ^ α ) and the application of their orthonormality relation, the calculation of the modified multipole moments inside the muffin-tin spheres is straightforward and results in q L α [ ρ α ] = ( 2 ℓ + 1 ) ! ! λ ℓ ∫ ​ 0 R α ρ L α ( r α )   i ℓ ( λ r α )   r α 2 d r α   . The summation of V λ I [ ρ ˘ α ] over all spheres α finally provides the contributions of the charges of all spheres to the interstitial potential. Equation 7 reduces to Eq. 40, since for the free-space Green function Eq. 39 the boundary value term disappears for r → ∞ .

 2.3 Interstitial Yukawa Potential

Suppose we had found a Fourier transformable pseudo-charge density ρ ¯ α inside the sphere consistent with the Yukawa potential produced outside the sphere, with coefficients ρ ¯ α ( K ) , and the Fourier series would converge rapidly throughout the periodic domain. Then we can find the Fourier coefficients of the pseudo-charge density, ρ ˜ ( K ) , by Eq. 8 and the solution of the modified Helmholtz Eq. 1 through Fourier transformation yields an algebraic equation from which we calculate the interstitial Yukawa potential, V λ I ( K ) = 4 π K 2 + λ 2 ρ ˜ ( K )   . 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 ρ ¯ α that replaces the true local charge density ρ ˘ α [see (11)] inside the muffin-tin sphere such that the Yukawa potential in the interstitial region, V λ I [ ρ ˘ α ] = V λ I [ ρ α ] − V λ I [ ρ I ] , due to the true charge density inside the sphere, is equal to the Yukawa potential in the interstitial region produced by the a priori unknown pseudo-charge density, V λ I [ ρ ¯ α ] = V λ I [ ρ ˘ α ] . From Eq. 41 we conclude that this is fulfilled if the modified multipole moments Eq. 42 of both charge densities, q L α [ ρ ˘ α ] = q L α [ ρ α ] − q L α [ ρ I ] and q L α [ ρ ¯ α ] , are equal. The modified multiple moments q L α [ ρ α ] are already known through Eq. 43. Next we determine the modified multiple moments q L α [ ρ I ] of the interstitial charge extended into the muffin-tin spheres and then construct the pseudo-charge density.

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 q L α [ ρ I ] in detail. We write ρ I relative to the sphere centre τ α and employ the Rayleigh expansion Eq. 30 to e i K ⋅ r α , which yields ρ I ( r ) = ∑ K ρ I ( K ) e i K ⋅ r α e i K ⋅ τ α = ∑ K ρ I ( K ) e i K ⋅ τ α ∑ L 4 π i ℓ j ℓ ( K r α ) Y L * ( K ^ ) Y L ( r ^ α )   . The modified multipole moments of ρ I in the sphere B R α ( τ α ) , defined analogously to Eq. 42, are q L α [ ρ I ] = ( 2 ℓ + 1 ) ! ! λ ℓ ∫ B R α ( 0 ) Y L * ( r ^ α ) i ℓ ( λ r α ) ρ I ( r α + τ α ) d r α = ( 2 ℓ + 1 ) ! ! λ ℓ 4 π i ℓ ∑ K ρ I ( K ) e i K ⋅ τ α Y L * ( K ^ ) ∫ ​ 0 R α i ℓ ( λ r α ) j ℓ ( K r α ) r α 2 d r α   . For K ≠ 0 , the latter integral becomes ∫ ​ 0 R α i ℓ ( λ r α ) j ℓ ( K r α ) r α 2 d r α = R α 2 K 2 + λ 2 ( K i ℓ ( λ R α ) j ℓ + 1 ( K R α ) + λ i ℓ + 1 ( λ R α ) j ℓ ( K R α ) ) = R α 2 K 2 + λ 2 ( λ i ℓ − 1 ( λ R α ) j ℓ ( K R α ) − K i ℓ ( λ R α ) j ℓ − 1 ( K R α ) )   . If K = 0 , observe that j ℓ ( 0 ) = δ ℓ 0 (Kronecker δ) and so the integral becomes δ ℓ 0 ∫ ​ 0 R α i 0 ( λ r α ) r α 2 d r α = R α 3 i 1 ( λ R α ) λ R α δ ℓ 0   . Both equations above can be derived by partial integration in two different ways and applying the identities in Ref. [18], d d r ( r − ℓ f ℓ ( r ) ) = ± r − ℓ f ℓ + 1 ( r ) for f ℓ = i ℓ with the plus sign and for f ℓ = j ℓ with the minus sign, and d d r ( r ℓ + 2 f ℓ + 1 ( r ) ) = r ℓ + 2 f ℓ ( r ) for both f ℓ = i ℓ and f ℓ = j ℓ . In conclusion, this yields q L α [ ρ I ] = δ ℓ 0 4 π R α 2 i 1 ( λ R α ) λ ρ I ( 0 ) + ∑ K ≠ 0 ( 2 ℓ + 1 ) ! ! λ ℓ 4 π i ℓ ρ I ( K ) e i K ⋅ τ α Y L * ( K ^ ) R α 2 λ 2 + K 2 ( K i ℓ ( λ R α ) j ℓ + 1 ( K R α ) + λ i ℓ + 1 ( λ R α ) j ℓ ( K R α ) )   .

Construction of Pseudo-Charge Density

We construct the pseudo-charge density by following the ansatz of Weinert, ρ ¯ α ( r α + τ α ) = ∑ L ρ ¯ L α ( r α ) Y L ( r ^ α ) = ∑ L Q L n α ( ∑ η = 0 n a η r α ν η ) Y L ( r ^ α )   , in which the radial dependence of the charge density is expressed in terms of a polynomial expansion up to degree ν n , which depends on atom α and angular degree ℓ , and otherwise use spherical harmonics for the angular part—this being the usual representation for charge densities in the muffin-tin region. As we will discuss in Section 2.3.3, it is beneficial to choose a η = ( − 1 ) n − η R α 2 ( n − η ) ( n η ) a n   for  η = 0 , … , n and ν η = ℓ + 2 η , where n is yet to be determined. As will become apparent later when we derive the coefficients Q L n α in Eq. 61, the coefficient a n cancels out in any relevant equation, like Eq. 54 or Eq. 55. With these choices of parameters and with the binomials theorem applied to ∑ η = 0 n ( − 1 ) n − η ( n η ) ( r α R α ) 2 η = ( ( r α R α ) 2 − 1 ) n it follows from the ansatz Eq. 51 ρ ¯ α ( r α + τ α ) = a n ( r α 2 − R α 2 ) n ∑ L Q L n α r α ℓ Y L ( r ^ α )   . The Fourier transform of this expression is then given by ρ ¯ α ( K ) = 1 | Ω | e − i K ⋅ τ α ∫ B R α ( 0 ) ρ ¯ α ( r α + τ α ) e − i K ⋅ r α d r α = 4 π | Ω | e − i K ⋅ τ α ∑ L ( − i ) ℓ Q L n α A ℓ n α ( K ) Y L ( K ^ )   , where A ℓ n α ( K ) = a n ∫ ​ 0 R α ( r α 2 − R α 2 ) n r α ℓ + 2 j ℓ ( K r α ) d r α , | Ω | is the volume of the periodic domain and we used the Rayleigh expansion Eq. 30 and the orthonormality relation of the spherical harmonics. With Proposition 1 in the appendix A ℓ n α ( K ) finally reduces to A ℓ n α ( K ) = a n ( − 2 ) n n ! R α ℓ + n + 2 j ℓ + n + 1 ( K R α ) K n + 1   . In the same way we derive the coefficients Q L n α : We insert Eq. 54 in the definition of the modified multipole moments Eq. 42 and use Eq. 1 for f ℓ = i ℓ and κ = λ to obtain q L α [ ρ ¯ α ] = ( 2 ℓ + 1 ) ! ! λ ℓ ∫ B R α ( 0 ) ρ ¯ α ( r α + τ α ) i ℓ ( λ r α ) Y L * ( r ^ α ) d r α = ( 2 ℓ + 1 ) ! ! λ ℓ Q L n α a n ∫ ​ 0 R α ( r α 2 − R α 2 ) n r α ℓ + 2 i ℓ ( λ r α ) d r α = ( 2 ℓ + 1 ) ! ! λ ℓ Q L n α a n ( − 2 ) n n ! R α ℓ + n + 2 i ℓ + n + 1 ( λ R α ) λ n + 1   , and thus, Q L n α = q L α [ ρ ¯ α ] λ ℓ + n + 1 a n ( − 2 ) n n ! R α ℓ + n + 2 ( 2 ℓ + 1 ) ! ! i ℓ + n + 1 ( λ R α )   . Since A ℓ n α ( K ) and Q L n α share the term a n ( − 2 ) n n ! R α ℓ + n + 2 , the term cancels out in the product Q L n α A ℓ n α ( K ) entering Eq. 55, Q L n α A ℓ n α ( K ) = j ℓ + n + 1 ( K R α ) K n + 1 ( 2 ℓ + 1 ) ! ! λ ℓ + n + 1 i ℓ + n + 1 ( λ R α ) q L α [ ρ ¯ α ]   , and setting ν = ℓ + n + 1 (in accordance with Section 2.3.3) this leads to ρ ¯ α ( K ) = 4 π | Ω | e − i K ⋅ τ α ∑ L ( − i ) ℓ j ν ( K R α ) K ν − ℓ ( 2 ℓ + 1 ) ! ! λ ν i ν ( λ R α ) q L α [ ρ ¯ α ]   Y L ( K ^ )   . For ρ ¯ α ( 0 ) we take the limit lim K → 0 j ν ( K R α ) K ν − ℓ = lim K → 0 K ℓ R α ν ( 2 ν + 1 ) ! !   , which is only different from 0 for ℓ = 0 and thus yields ρ ¯ α ( 0 ) = 4 π | Ω | ( λ R α ) ν ( 2 ν + 1 ) ! ! i ν ( λ R α ) q 00 α [ ρ ¯ α ]   . 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 ρ ˘ α , q L α [ ρ ˘ α ] = q L α [ ρ α ] − q L α [ ρ I ] calculated from the modified multipole moments q L α [ ρ α ] Eq. 43 and q L α [ ρ I ] Eq. 50 of ρ α and ρ I respectively, in the sphere B R α ( τ α ) .

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

In Section 2.3.2, we have set a η by Eq. 52, ν η = ℓ + 2 η and ν = ℓ + n + 1 without having determined n yet. Here we motivate our choices and finally determine a proper n.

With our choices for a η and ν η we have eradicated the sum in ansatz Eq. 51 and derived the much simpler form Eq. 54. The function r α ↦ ( r α 2 − R α 2 ) n itself and all its first n − 1 derivatives with respect to r α are equal to zero at r α = R α . Consequently, we ensure smoothness on the boundary of the sphere, d k d r α k ρ ¯ L α ( r α ) | r α = R α = 0   ∀ k = 0 , … , n − 1 . Note that this includes localization of ρ ¯ α in B R α ( τ α ) and thus the pseudo-charge density equals the charge density in I (condition 2 in Section 2.1).

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 n − 1 -fold differentiable function in combination with the differentiation rule for a Fourier transform, one can show that the Fourier coefficients for K → ∞ go faster to zero than 1 K n − 1 , i.e., we obtain the fastest convergence of the Fourier series for large 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 Q L n α A ℓ n α ( K ) / q L α [ ρ ¯ α ] , where in the Coulomb case Q L n α A ℓ n α ( K ) = j ℓ + n + 1 ( K R α ) K n + 1 ( 2 ℓ + 1 ) ! ! ( 2 ℓ + 2 n + 3 ) ! ! R α ℓ + n + 1 q L α [ ρ ¯ α ]   . His reasoning is based on the zeros of the K-dependent function Q L n α A ℓ n α ( K ) / q L α ( ρ ¯ α ) . Viewed as a function of K, however, our factor Q L n α A ℓ n α ( K ) / q L α [ ρ ¯ α ] differs from Weinert’s one only by a multiplicative constant. Thus Weinert’s arguments apply here as well. With Proposition 2 in the Appendix it follows that this multiplicative constant is smaller than 1 for λ > 0 . We see this confirmed in Figure 1, which shows the K R α -dependence of the factor Q L n α A ℓ n α ( K ) / q L α [ ρ ¯ α ] in the Yukawa Eq. 62 and Coulomb Eq. 68 cases for several combinations of n and angular degree ℓ — it reveals a smaller amplitude in the Yukawa case. In agreement with Ref. [12], in Figure 1 we make two observations: (i) Q L n α A ℓ n α / q L α [ ρ ¯ α ] as a function of K R α has larger oscillations for smaller n, and (ii) for fixed n, the largest contribution to the Fourier series comes from K R α less than the first zero of the Bessel function j ℓ + n + 1 . Since we deal with a finite number of 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 K max , the factor Q L n α A ℓ n α / q L α [ ρ ¯ α ] must be small for K > K max , which limits ℓ + n + 1 to a certain value, since the first zero is pushed towards infinity for growing ℓ + n + 1 , as can be seen from a comparison between the pink and yellow, or the blue and purple lines in Figure 1. From this arises Weinert’s criterion for choosing n ( ℓ ) , which we adopt here:

• Choose ν ∈ ℕ such that the first zero of j ν ( z ) is approximately equal to ( K R α ) max .

• Then n ( ℓ ) is fixed by the relation ν = ℓ + n + 1 .

The factor Q L n α A ℓ n α / q L α [ ρ ¯ α ] for several ( ℓ , n ) in the Yukawa and Coulomb cases with R α = 1 and λ = 2 .

Note that in this method ℓ is compensated by n in such a way that ν is de facto not depending on ℓ . Since the discretization of 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. Bulk-Case Yukawa Potential

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

Output: Yukawa potential V λ solving the modified Helmholtz Eq. 1 with periodic boundary conditions.

 Pseudo-Charge Density ρ ˜ ← ρ

1: Modified multipole moments q L α [ ρ I ] of the interstitial charge density in B R α ( τ α ) . ▷ Eq. 50

2: Modified multipole moments q L α [ ρ α ] of the muffin-tin charge density in B R α ( τ α ) . ▷ Eq. 43

3: Modified Multipole Moments q L α [ ρ ¯ α ] = q L α [ ρ α ] − q L α [ ρ I ] of ρ ¯ α

4: Sphere-localized part ρ ¯ α ( K ) of the pseudo-charge density. ▷ Eqs. 63 and 65

5: Pseudo-charge density ρ ˜ ( K ) . ▷ Eq. 8 Interstitial potential V λ I ← ρ ˜

6: Interstitial potential V λ I ( K ) . ▷ Eq. 44 Muffin-tin potential V λ MT ← ρ MT , V λ I

7: Boundary terms V L I ( R α ; τ α ) of muffin-tin potential. ▷ Eq. 31

8: Radial parts V L α ( r α ) of muffin-tin potential. ▷ Eq. 38

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, q L α [ ρ I ] Eq. 50 and q L α [ ρ α ] Eq. 43 respectively, and in the Fourier components of the pseudo-charge density ρ ¯ α ( K ) Eqs. 63 and 65. The integer ν 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 4 π / K 2 + λ 2 ) that substitutes 4 π / K 2 . Since the 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 r ℓ and 1 / r ℓ + 1 in these quantities are substituted by the modified spherical Bessel functions i ℓ ( λ r ) and k ℓ ( λ r ) , respectively, and the prefactor 4 π / ( 2 ℓ + 1 is substituted by 4 π λ . The interstitial potential on the boundary of the spheres, V L I ( R α ; τ α ) Eq. 31, only changes indirectly through the changed values of V λ I ( K ) .

 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 1 / λ . Like Weinert’s pseudo-charge method, our extension is based on the concept of the non-uniqueness of the multipole expansion as well as on the Dirichlet boundary value problem applied to a sphere. The difference between the modified Helmholtz equation and the Poisson equation lies solely in the radial behavior and thus the homogeneous solutions to the radial part of the differential equation are modified spherical Bessel functions instead of polynomial functions. The consequence is a different radial behavior of the Green function, resulting in the screening of the potential, which is now expanded in modified multipole moments, and this in turn affects the pseudo-charge density. Furthermore, the modification of the multipole moments implies that the modified monopole is not connected anymore to the total charge. 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 and thus we can therefore best choose the same integer parameters for convergence. Finally we layed out the minor changes necessary to change an implemented method for solving the Poisson equation available to an implementation for solving the modified Helmholtz equation.

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.

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].

 References Yukawa H . On the interaction of elementary particles. I. Proc Physico-Math Soc Jpn (1935) 33:48–57. 10.11429/ppmsj1919.17.0_48 Hückel E Debye P . The theory of electrolytes. I. lowering of freezing point and related phenomena. Phys Z (1923) 24:185–206. Thomas L . The calculation of atomic fields. Math Proc Cambridge Philos Soc (1927) 23:542–8. 10.1017/S0305004100011683 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. 10.1007/BF01351576 Hohenberg P Kohn W . Inhomogeneous electron gas. Phys Rev (1964) 136:B864–71. 10.1103/PhysRev.136.B864 Kohn W Sham L . Density function theory. J Phys Rev (1965) 140:A1133–8. 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. 10.1103/PhysRevB.83.235118 Massidda S Posternak M Baldereschi A . Hartree-Fock lapw approach to the electronic properties of periodic systems. Phys Rev B (1993) 48:5058–68. 10.1103/PhysRevB.48.5058 Kerker G . Efficient iteration scheme for self-consistent pseudopotential calculations. Phys Rev B (1981) 23:3082. 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. 10.1103/PhysRevB.102.195138 Kim J Gulans A Draxl C . Robust mixing in self-consistent linearized augmented planewave calculations. Electron Struct (2020) 2:037001. 10.1088/2516-1075/ababde Weinert M . Solution of Poisson’s equation: beyond Ewald-type methods. J Math Phys (1981) 22:2433–9. 10.1063/1.524800 Eyert V . The plane-wave based full-potential ASW method. Berlin, Heidelberg, Germany: Springer (2013). p. 113–72. 10.1007/978-3-642-25864-0˙4 Drittler B . KKR-Greensche Funktionsmethode für das volle Zellpotential. [PhD thesis]. Aachen (Germany): RWTH Aachen (1991). 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: O2 molecule. Phys Rev B (1981) 24:864–75. 10.1103/PhysRevB.24.864 FLEUR. The Jülich FLAPW code family. Available from: http://www.flapw.de/ (2020). 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). Arfken G Weber H Harris F . Mathematical methods for physicists: a comprehensive guide. Amsterdam, Netherlands: Elsevier Science (2013). Jackson JD . Classical electrodynamics. 3rd ed. New York, NY: Wiley (1999). Jülich Supercomputing Centre. JURECA: modular supercomputer at Jülich supercomputing centre. J Large-Scale Res Facil (2018) 4:33. 10.17815/jlsrf-4-121-1 Hinzen M Napoli ED Wortmann D Blügel S . Solution to the modified helmholtz equation for arbitrary periodic charge densities (2020). APPENDIX

Proposition 1. Let R > 0 , κ > 0 and f ℓ be either the spherical Bessel function of the first kind, f ℓ = j ℓ , or the modified spherical Bessel function of the first kind, f ℓ = i ℓ . Then ∫ 0 R ( r 2 − R 2 ) n r ℓ + 2 f ℓ ( κ r ) d r = ( − 2 ) n n ! R ℓ + n + 2   f ℓ + n + 1 ( κ R ) κ n + 1 holds for all n ∈ ℕ 0 and all ℓ ∈ ℕ 0 .

Proof by mathematical induction on n. The base case n = 0 follows immediately from the differentiation property Eq. 49 of the functions j ℓ and i ℓ , by ∫ 0 R r ℓ + 2 f ℓ ( κ r ) d r = κ − ℓ − 3 ∫ 0 κ R r ℓ + 2 f ℓ ( r ) d r = κ − 1 R ℓ + 2 f ℓ + 1 ( κ R ) . The induction step uses the property for n − 1 and ℓ + 1 and partial integration with u ( r ) = ( r 2 − R 2 ) n and v ′ ( r ) = r ℓ + 2 f ℓ ( κ r ) , to derive the statement for n and ℓ . Let F ( n , ℓ ) = ∫ ​ 0 R ( r 2 − R 2 ) n r ℓ + 2 f ℓ ( κ r ) d r be the left-hand side of Eq. 1. Then F ( n , ℓ ) = [ ( r 2 − R 2 ) n κ − 1 r ℓ + 2 f ℓ + 1 ( κ r ) ] 0 R − ∫ 0 R n ( r 2 − R 2 ) n − 1 2 r κ − 1 r ℓ + 2 f ℓ + 1 ( κ r ) d r       = − 2 n κ − 1 F ( n − 1 , ℓ + 1 ) = ( − 2 ) n n ! R ℓ + n + 2   f ℓ + n + 1 ( κ R ) κ n + 1 , where we used the induction hypothesis in the last equation, and thus the proposition follows. ∎

Proposition 2. For nonnegative λ and R, and ν ∈ ℕ i ν ( λ R ) λ ν ≥ R ν ( 2 ν + 1 ) ! ! holds, with equality if and only if λ = 0 or R = 0 .

 Proof. Since i ν ( x ) = i − v j ν ( i x ) and j ν has the expansion [18] j ν ( x ) = x ν ( 2 ν + 1 ) ! ! ∑ s = 0 ∞ ( − 1 ) s s ! ( ν + 3 2 ) s ( x 2 ) 2 s   , where ( ⋅ ) s is the Pochhammer symbol defined for general a ∈ ℝ and s ∈ ℕ 0 by ( a ) 0 = 1   ,   ( a ) s = a ( a + 1 ) ⋯ ( a + s − 1 )   , i ν can be expanded by i ν ( x ) = x ν ( 2 ν + 1 ) ! ! ∑ s = 0 ∞ 1 s ! ( ν + 3 2 ) s ( x 2 ) 2 s   . i ν ( λ R ) / λ ν thus becomes i ν ( λ R ) λ ν = R ν ( 2 ν + 1 ) ! ! ∑ s = 0 ∞ 1 s ! ( ν + 3 2 ) s ( λ R 2 ) 2 s and it remains to show that ∑ s = 0 ∞ 1 s ! ( ν + 3 2 ) s ( λ R 2 ) 2 s ≥ 1   . The first term of the sum for s = 0 is always 1 irrespective of the values of λ and R. The other summed terms for s > 0 are positive, if both λ and R are positive, and zero, if one of the two variables is zero. ∎