Electrostatic Interaction of Point Charges in Three-Layer Structures: The Classical Model

: Electrostatic interaction energy W between two point charges in a three-layer plane system was calculated on the basis of the Green’s function method in the classical model of constant dielectric permittivities for all media involved. A regular method for the calculation of W ( Z , Z (cid:48) , R ) , where Z and Z (cid:48) are the charge coordinates normal to the interfaces, and R the lateral (along the interfaces) distance between the charges, was proposed. The method consists in substituting the evaluation of integrals of rapidly oscillating functions over the semi-inﬁnite interval by constructing an analytical series of inverse radical functions to a required accuracy. Simple ﬁnite-term analytical approximations of the dependence W ( Z , Z (cid:48) , R ) were proposed. Two especially important particular cases of charge conﬁgurations were analyzed in more detail: (i) both charges are in the same medium and Z = Z (cid:48) ; and (ii) the charges are located at different interfaces across the slab. It was demonstrated that the W dependence on the charge–charge distance S = (cid:113) R 2 + ( Z − Z (cid:48) ) 2 differs from the classical Coulombic one W ∼ S − 1 . This phenomenon occurs due to the appearance of polarization charges at both interfaces, which ascribes a many-body character to the problem from the outset. The results obtained testify, in particular, that the electron–hole interaction in heterostructures leading to the exciton formation is different in the intra-slab and across-slab charge conﬁgurations, which is usually overlooked in speciﬁc calculations related to the subject concerned. Our consideration clearly demonstrates the origin, the character, and the consequences of the actual difference. The often used Rytova–Keldysh approximation was analyzed. The cause of its relative success was explained, and the applicability limits were determined.


Introduction
Electrostatic interactions in layered systems have been extensively studied in the literature, making allowance, if necessary, for the spatial and temporal dispersion of the dielectric functions ε(k, ω) in the constituent media . Here, k and ω are the transferred wave number and frequency, respectively. Theoretical papers devoted to the problem can be grouped according to: (i) the number of layers in the system (two, three, or the infinite number); or (ii) which interaction is considered. The latter means either the test-charge (in rare cases, test-dipole) electrostatic self-energy arising due to the existence of interface(s) between the media with different dielectric permittivities, which is usually called the image force energy [25][26][27][28][29][30][31][32][33], or the interaction energy between the charges (or dipoles) located in the same layer or in different layers [6,8,[34][35][36][37][38][39][40][41][42][43].
In the papers cited above, it was demonstrated that the effects of spatial and temporal dispersion are important and sometimes crucial when describing various phenomena governed by the layered character of the system. Nevertheless, certain phenomena in the layered structures, e.g., exciton-related ones, can be understood even in the ε = const (constant-dielectric-permittivity, CDP) approximation. The latter brings about some interesting results that are not eliminated by the account of the spatial dispersion effects (the screening non-locality) [11,13,15,22,[44][45][46][47][48][49][50][51].
Note that the image forces for charges in three-layer structures were studied rather well theoretically, whereas our knowledge of charge-charge interactions in those objects remains unsatisfactory, although the latter problem is not less important for applications [35,[52][53][54][55][56][57][58]. For instance, the interaction between ions adsorbed on a solid surface or electrons located above a liquid helium surface is usually (and, in some cases, quite correctly!) assumed to be a classical Coulombic one. As for electrons above helium films, the distortion of the Coulomb charge-charge interaction law (up to its conversion into the dipole-dipole one) due to the existence of a substrate with its own dielectric permittivity that is much higher than that of helium was recognized long ago [50]. Nevertheless, a proper comprehensive analysis of the phenomenon concerned is still lacking.
In this article, making use of calculation schemes used earlier [6,11,15], we obtain explicit formulas for the interaction W between two charges located in a three-layer system composed of two semi-infinite media separated by a flat interlayer, with all three media having different electrostatic properties (dielectric permittivities). The results obtained demonstrate that, in any layered system, the polarization of the media violates the Coulombic dependence of W on the distance S between the charges, i.e., the quantity W is not inversely proportional to S! Earlier we drew the same qualitative conclusion for the simplest inhomogeneous system consisting of two different semi-infinite homogeneous media separated by a plane interface [59]. However, in that case, the application of textbook electrostatic formulas was enough. The three-layer configuration is much more complicated from the mathematical viewpoint. In this paper, we show how the interaction energy W and its long-range asymptotics can be accurately calculated with a required accuracy in the framework of the CDP model and propose rather useful analytical approximations. To make the paper less cumbersome, we did not consider the cases when a proper account of the spatial dispersion is necessary.

Formulation of the Problem and Preliminary Remarks
Let two point-like test charges Q and Q be located in a three-layer system composed of two semi-infinite media 1 and 3 (hereafter the covers) separated by flat interlayer 2 (hereafter the slab) of the thickness L. All three media are considered to be spatially homogeneous and characterized by their own bulk dielectric functions ε i (k, ω), i = 1, 2, 3. We have to determine the electrostatic interaction energy between those charges and its dependence on various problem parameters. In Appendix A, the general solution of the problem in terms of Green's functions is presented. However, no analytical results were obtained in the general case (see [6,11,15,43] and the references therein), whereas the computational efforts to obtain illustrative and useful results were substantial. Our current treatment is carried out in the framework of a much simpler classical electrostatic model (the CDP one), when the spatial dispersion of the dielectric functions ε(k, ω)s is neglected, i.e., ε i (k, ω) = ε i (0, ω) for each ith medium. Since only static phenomena (fixed charge locations) are dealt with, the energy W dependence on ω becomes irrelevant, so that all involved dielectric functions are hereafter treated as constants, and the ω-dependence is omitted from consideration. This approach allowed us to find an effective method to solve the problem, which we believe to be also applicable in more general cases.
In what follows, we have to specify the coordinate frame. There are two evident choices. Each of them has its own specific advantages. In one of them, the charge coordinates Z and Z are reckoned from the interface between media 1 and 2 perpendicularly to both interfaces, with the Z-values being negative in medium 1 (Z < 0) and positive in media 2 (0 < Z < L) and 3 (Z > L). The longitudinal (along the interfaces) projection of the segment between the charges equals R, so that the actual distance is This choice is convenient if we are mainly interested in the processes running in either of the covers and intend to study their dependence on the interlayer thickness L. For the sake of completeness and illustration, the most general solution of the problem given in Appendix A is presented just in this reference frame.
The other choice consists in identifying the middle plane of interlayer 2 with the XY plane. Then, e.g., Z < −L/2 in medium 1, −L/2 < Z < L/2 in medium 2, and Z > L/2 in medium 3. Equation (1) for the distance between the charges remains the same in this reference frame. The advantage of this choice consists in that certain symmetry properties of the problem formulation and the results (e.g., while substituting ε 1 ↔ ε 3 ) are traced more easily. Because of the latter reason, this reference frame was selected for the further analysis (see Figure 1). Below, to make the equations more compact, all length-related problem parameters-these are R, S, Z, and Z -are almost everywhere normalized by the characteristic length of the problem, the interlayer width L. Hereafter, the normalized quantities will be designated by the same letters, but in the lower case: r, s, z, and z . Then, the general expression for W QQ can be factorized in the form

Zʹ
which we use as a definition of the dimensionless charge-charge interaction energy w. Strictly speaking, the width L has to be substantially larger than the atomic scale, which is not always the case in heterostructures. However, even if this requirement is not satisfied, the formally macroscopic dielectric approach can be used, at least as a guide. It is so because it is known that many properties, including the superconducting phase transition or the screening ability, do not disappear in the limit of atomically thin solid films. Hence, for qualitative purposes, the interlayers may be treated as continuous two-dimensional solids [94][95][96][97][98]. For instance, the macroscopic treatment of the dielectric function is often applied even to graphene, which is a one-layer material [99,100]. Of course, in the latter case, the graphene "dielectric constant" is only an effective screening parameter of a two-dimensional sheet [101] and should not coincide with the bulk graphite value. However, the parameter concerned together with the dielectric constants of the environment may be used in the classical description of layered structures. Surely, the account of the dielectric function spatial dispersion can substantially improve the treatment.
All analytical calculations were carried out in the Maple environment.

General Formulas for w(z, z ) in the Dispersionless Approximation
In the adopted normalized reference frame, the z-coordinates in the slab are confined within the interval − 1 2 , 1 2 , whereas the semi-space z < − 1 2 is occupied by medium 1, and the semi-space z > 1 2 by medium 3 (see Figure 1). After performing all required normalizations, one can find that the electrostatic interaction energy between two point-like charges with the arbitrary coordinates z and z and separated by the lateral (along the interfaces) distance r can be written in the general form (cf. Equation (A1)) being a consequence of the transition from the three-dimensional Fourier transform to the Hankel transform in the actual case of the circular symmetry of the problem [102,103] (all media are assumed to be isotropic).
Here, q is a dimensionless integration variable, J 0 (x) is the Bessel function of the zeroth order, and the notation ε implies the whole set of dielectric constants {ε 1 , ε 2 , ε 3 }. An additional bonus of introducing the factor Sub(q, z, z ; ε) = qD(q, z, z ; ε) is the fact that we get rid of the inherent q −1 -singularity of Green's function D(q, z, z , ε) in the two-dimensional momentum q-space [35,104,105] (see Appendix A). Equation (3) is remarkable per se, because only the Bessel function in the integrand depends on the lateral distance r between the charges, irrespective of their locations in the system. At the same time, the introduced factor Sub(q, z, z ; ε), which is referred to as subintegrand throughout this paper, depends on all other problem parameters, but r. Therefore, according to the properties of the Hankel transform, we obtain a paradoxical at first glance situation when the function Sub(q, z, z ; ε), being independent of r, unambiguously determines the dependence of w on this variable.
In turn, Sub(q, z, z ; ε) can be represented as the ratio The denominator includes the parameter and describes the interference of electric fields induced by polarization charges located at both interfaces. It is examined in Appendix B in more detail. The denominator is the same for any arrangement of charges in the heterostructure. The specificity of charge locations is reflected in the numerator Num(q, z, z ; ε). This quantity can be presented in the form where the non-negative-valued functions f i (z, z ), i.e., f i (z, z ) ≥ 0, can acquire the following values: One can see that these values can be paired with respect to the specific combinations ± z ± z with the pair components differing by 2. The number of terms, N, in the sum in Equation (7) and the specific combinations involved (Equation (8)) are determined by the arrangement of charges. In Appendix C, we present the full sets of summands {F i } in Equation (7) that are required to calculate the interaction energy between two point-like charges with arbitrary locations in the three-layer system concerned.
For certain z and z , e.g., z = z , some of the involved functions f i (z, z ) can coincide, and the number N becomes smaller. The calculation of the integral in Equation (3) is a trivial task if the parameter θ = 0 (in this case Den(q, ε) = 1), and it can be easily carried out analytically (see Appendix D). According to Equation (6), this case corresponds to either ε 2 = ε 1 or ε 2 = ε 3 , so that the three-layer problem reduces to the two-layer one [59,106]. If ε 1 = ε 2 = ε 3 , we obtain a homogeneous medium over the whole space. However, if θ = 0, the presence of the Bessel function, which can hardly be substituted by anything more simple within the interval [0, ∞) of its argument, makes the integration much more difficult. All attempts to analytically approximate the integral in Equation (3) are reduced to nothing else but various approximations of the function Num(q, z, z , ε) in Equation (4).

Regular Method of Integral Calculation
In the general case θ = 0, the calculation procedure becomes more complicated owing to the second term in Equation (5). Nevertheless, let us take into account that both |θ| < 1 [we exclude the "weird" cases (ε 1 , ε 2 , ε 3 ) → ∞ from consideration] and exp (−2q) < 1 at q > 0, so that the absolute value of this term is practically always less than unity, and the denominator itself never equals zero. All summands in the function F(q, z, z ; ε) are finite. Furthermore, most of them tend to zero at q → ∞. Therefore, we conclude that integral (3) converges for all r > 0 and for all relevant ε i .
Let us consider the second term in Den(q, ε) as a small parameter and expand Den −1 (q, ε) into a convergent series Then, taking into account Equation (7), we obtain an infinite series of terms, each of which, according to Equation (A35), is integrable analytically. As a result, the energy w(z, z , r; ε) is expressed as a finite sum of several infinite series If θ > 0, Equation (11) represents the alternating series with monotonically decreasing terms. If θ < 0, the series terms can be majorized by |θ| j / (2j). Accordingly, we can easily evaluate the corresponding absolute error for each series. In both cases, the sums are convergent, but the convergence rate is much faster for θ > 0.
Nevertheless, it is not worthwhile to evaluate each series separately up to a certain relative accuracy and thereafter sum up the results, because the summands may be of different signs and strongly compensate one another. Hence, we used the following algorithm. Each next term in the total sum in Equation (10) was selected as having a maximum magnitude among the remaining terms (recall that they monotonically decrease in each series) in all series. The absolute error of this procedure was determined by summing up the absolute errors of all series.
We would like to draw attention to the following points. Firstly, the described expansion into a series can be made not simultaneously, but step by step. Let us examine the ith term in Equation (7), i.e., the integral Multiplying its integrand by 1 = 1 + θ exp (−2q) − θ exp (−2q) and using Equation (A35), this expression can be transformed into another one The first term does not need any additional simplification and can be placed to what is referred to as the Coulomb tail. This tail will accumulate Coulomb-like terms produced by every summand in Equation (7). The second term in Equation (13) can be further iteratively processed using the same procedure, so that new terms will appear to be incorporated into the Coulomb tail. This trick can be applied to reduce the number of terms in the sum in Equation (7), for which the difference between the arguments of exponential functions is a multiple of 2q (cf., e.g., the terms F 1 and F 4 , and F 2 and F 3 for w 11 in Appendix C) and simplify the following summation procedure. Equation (8) testifies that such a situation can really take place. Below, we refer to the indicated procedure as term grouping. For specific zand z -values, a group may include several F i terms from the sum in Equation (7). Their power exponents belong to the set in Equation (8) and differ by multiples of 2. The total contribution of this group to w(z, z , r; ε) will include a few terms similar to the first summand in Equation (13) and included into the Coulomb tail, and a single common integral such as the second summand in the same equation. Sometimes, as shown below, the common integral may vanish altogether, which significantly simplifies the calculation.
Secondly, the described procedure of series summation is governed by the required accuracy of the final result. Therefore, the number of calculated terms incorporated into the Coulomb tail is not known in advance. It depends on the problem parameters and in certain cases can be rather large. Evidently, any analytical result in a compact form is out of the question in this case. However, one may talk about finite analytical expressions that would approximate the integral in Equation (3) rather well. In Equation (13), the first term is an exact analytical contribution to the calculated result, and the second one can be regarded as a "correction". The argument of the exponential function in the integrand of the second term has an increment of 2 in comparison with initial Equation (12), so that the convergence of this "correction" is better. Hence, a more accurate approximation can be obtained for this quantity. Below, we demonstrate that the trick of extracting "exact" contributions to the final result, i.e., the Coulomb-tail term may be quite reasonable. However, the procedure can be repeated again and again, and the decision when to terminate it is a rule of thumb.
Nevertheless, some considerations of a general character can be made even at this early stage. By analyzing the formulas in Appendix C, one can see that each set of summands contains a term with f (z, z ) = |z − z | (we put it in the first place, F 1 ). This term describes the direct Coulomb interaction between the charges, but modified (screened) by the interface polarization charges. One can verify (we demonstrate it below in the case z = z ) that, if the charges approach infinitesimally close to each other (if the problem specification allows it), the bare (unscreened) Coulomb contribution will dominate in their interaction. In other words, the charges cease to feel the three-layer character of the system. The interaction between them will be Coulombic, w(s → 0) ∼ s −1 , with the coefficient of proportionality being equal to the coefficient in F 1 . On the other hand, one can easily obtain the following remarkable result. If r → ∞ provided that z and z are fixed (i.e., s → ∞), then In the second row, we made the substitution rq → q in the integral. Now, it is easy to verify that for any {F i } set. As a result, we obtain that lim s→∞ w(z, z , r; ε) = 2 irrespective in which media the charges are located, so that any information about the (existing!) interlayer disappears. The same formula describes the interaction between two point-like charges located at the single planar interface between homogeneous media in a two-layer system, where the interlayer is absent [59]. A similar character of w-asymptotics (w ∼ s −1 , although with different coefficients) at s → 0 and s → ∞ testifies that the product sw is an excellent quantity to be analyzed while studying the behavior of the interaction energy w itself. It can be considered as the coefficient C in the representation of w in the Coulomb-like form w = C s .
However, the variety of problem parameters makes the general analysis too cumbersome. That is why hereafter we confine ourselves to two representative charge arrangements in the system concerned: either the charges are located in the same plane parallel to the interfaces (i.e., z = z ), or they are positioned at the opposite sides of the slab (i.e., z = − 1 2 and z = 1 2 ). Furthermore, those configurations are the most interesting from the practical viewpoint.
For numerical illustration, we selected toy three-layer heterostructures with ε-values quoted in Table 1. In the next section, we show that, besides the "long-range" characteristic value (Equation (15)), the parameter C in Equation (17) takes on the following "short-range" characteristic values 1 ε 1 , 2 For three non-symmetric (ε 1 = ε 2 = ε 3 ) structures with comparable ε i s (A-C), we chose layer dielectric constants equal to 1, 2, and 5 taken in various combinations. Those values were chosen to ensure that any value among Equations (15) and (18) does not coincide with another one. A value of 1.5 for symmetric (ε 1 = ε 3 = ε and ε 2 = ε ) structures D and E was selected, because this ratio is typical of currently studied heterostructures [62,72,[107][108][109][110]. Large values of 10 and 100 were taken for symmetric structures F to I following the seminal work by Keldysh [48], who discussed sandwiches containing polar solids with large dielectric constants. All ε sets mentioned in Table 1 are also indicated in Figure 2.
The lettered points correspond to the parameters of illustrative toy structures (see Table 1). Symmetric structures lie on the dashed diagonal.

The Particular Case z = z . Series Calculations
In this and several following sections, we consider the case, when the both charges are located at the same vertical level (z = z ≡ z) in the three-layer system. Then, the normalized distance s between them (see Equation (56)) becomes equal to the parameter r, which describes their relative interval in the lateral plane (see Figure 1). Therefore, a quantity to be analyzed is the product rw. Now, only three {F i } sets are relevant for the calculations: • for w 11 (z ≤ − 1 2 ): • for w 33 ( 1 2 ≤ z): In Figure 3, the calculation results obtained for the chosen three non-symmetric toy configurations are displayed. The presented pattern corresponds to a pair of charges separated by the spacing r and moving perpendicularly to the interfaces. The figure demonstrates that the interaction energy between the charges, if the distance r = 0, remains a continuous and smooth function, even when the charges cross either of the interfaces. However, this smoothness is a specific feature of the charge arrangement z = z . Otherwise, the polarization charges induced at the interface by a charge located near the interface would result in a cusp in the dependence w(r, z = z ) when the other charge crosses the interface. However, if z = z , there are no polarization charges at the interface at the moment when the both charges cross it simultaneously(!): the problem symmetry dictates that the electric field vectors at the interface must lie in the horizontal plane at this moment. Nevertheless, small remnants (or, better to say, precursors) of the indicated non-smoothness can be observed at small zs, i.e., when the charges are rather close to each other.  Strictly speaking, the realistic motion of charges perpendicularly to the interfaces should be considered making allowance for the change of the system's total electrostatic energy. The latter, besides the energy of interaction between the charges, also includes the image force energy (self-energy). However, the charge motion across the interfaces presented here is virtual. Namely, below, we assume that the charge positions with respect to the interfaces is fixed (z = const, z = const), so that we are interested in the r-dependence of the charge-charge interaction energy w(r). In this case, the variation of the image force energy does not matter, and the change in the total electrostatic energy is determined by that of w(r).
In Section 3, it is shown how the functions w(r) and Sub(q) are related by the Hankel transform (see Equation (3)). Figure 4 illustrates the corresponding pairs of functions rw(r) and Sub(q) in a wide interval of the charge coordinate z in toy structure A (this is a counterpart of Figure 3a). For this {ε i } combination, the depicted dependences do not overlap, but transform monotonically as z varies, which allows all typical features to be traced more or less clearly.
We draw attention to a close resemblance of the plots in Figure 4a,b. This similarity is true for the central section of the plots, being very useful, because it allows one to predict the behavior of the rw(r) (and, hence, w(r)) dependence at the qualitative level by considering Sub(q) without integration. Vice versa, one can judge about the behavior of Sub(q) on the basis of the dependence rw(r). This property dictated by the Hankel transform relation [102,103] between w(r) and D(q) (see Equation (3)) will be used in the further analysis. However, as follows from Figure 4c,d, the functional dependences of Sub(q) on q at q → 0 and rw(r) on 1/r at r → ∞ are different. Since rw(r) rather than Sub(q) is a matter of major concern, let us consider it in more detail. In the intervals 0.4 < |z| < 0.6, z is varied with a step of 0.01 in order to illustrate the transformation of curves when crossing the interfaces. (c) Short-wave (q → 0) asymptotics of the dependence Sub(q) for toy structure A and various zs from −10 to 10. The q-axis scale is linear. (d) Long-range (r → ∞) asymptotics of the dependence rw(r) for toy structure A and various zs from −10 to 10. The r-axis scale is reciprocal. The quantities q, z, z , and r are dimensionless.
The main features of the rw(r) dependence are as follows.

•
One can see that the quantity rw (in essence, the strength of the normalized Coulomb interaction in a three-layer system) strongly depends on z and r. This fact reflects the many-body character of the image forces [111,112] crucially distorting the parent electrostatic interaction in the inhomogeneous (in particular, layered) media.

•
For infinitesimally close charges located in the covers' bulk, i.e., not near the interfaces, the product rw approaches the appropriate values ε −1 i that characterize the corresponding bulk media i. If the charges are at the interface between media i and j, this value equals 2/ ε i + ε j , as in much simpler two-layer structures [59]. Thus, when two infinitesimally close charges (r → 0) cross an interface (z = ± 1 2 ) between media i and j, the value of the product rw changes in a step-like manner. To be more precise, it demonstrates two jumps: from ε −1 i to 2/ ε i + ε j and further to ε −1 j . (Note that such an analysis would be impossible for the "bare" (truly Coulombic) function w(r) possessing a conventional pole singularity at r → 0.) At r = 0, the product rw-and, hence, the function w(r) itself-changes smoothly [as was illustrated in Figure 3a. This transformation was considered in more detail in [59] for a two-layer system. The analysis remains valid for the three-layer system, because two infinitesimally close charges near either of the interfaces "cease to feel" the other interface (as was in the case of the charges in the bulk, where the specific Coulombic singularity "masks" all other interactions in the system), and the three-layer problem transforms into the two-layer one. In the vicinity of interfaces (z = ± 1 2 ), we plotted the dependences rw(r) corresponding to various zs with a small step in order to illustrate how the smoothness of the dependence rw(r = 0, z) matches the jumps in the dependence rw(r = 0, z).

•
At r → ∞, all dependences tend to the same limit 2/ (ε 1 + ε 3 ) (see Equation (16)). In terms of dimensional variables, the limit r → ∞ corresponds to both R → ∞ at a fixed L and L → 0 at a fixed R. Since Z = zL, the latter interpretation can be regarded as a consequence of two processes: (i) the disappearance of the slab (we leave the ideal conductor case ε 2 = ∞ beyond the consideration), so that the three-layer system transforms into a two-layer one; and (ii) the deposition of both charges onto a single interface between two media (1 and 3). From the qualitative viewpoint, this limiting case means that there are no electric force lines in the slab.

•
Even for charges located in the depth of any medium, their interaction deviates from the "true" Coulombic dependence (∼ r −1 ) already at rather short distances from each other, i.e., the influence of other media matters.
Some of the rw-dependences are non-monotonic. However, we emphasize (see Equation (17)) that the quantity rw is only a coefficient in front of r −1 . The "bare" interaction energy w remains to be a monotonically decreasing function of r, because it has to satisfy the Laplace equation where ∆ is the Laplace operator, in the bulk of each medium, since the polarization charges in heterostructures arise only at the interfaces.
The corresponding dependences for toy structures B and C have the same characteristic features. However, the curves in the relevant plots strongly intersect one another and can lead to a confusion. Therefore, those sets are not presented. However, the picture becomes monotonic again and convenient for analysis in the case of symmetric structures if we consider either of the symmetrical semi-spaces z ≥ 0 or z ≤ 0. Figures 5 and 6 illustrate the corresponding Sub-rw pairs, as well as their short-wave and long-range asymptotics, for various zs in toy structures D and E, respectively. One can see that the similarity between the curves in the pairs is not so clear, but it still exists. The main difference between the curves displayed in Figures 5 and 6 and those depicted in Figure 4 consists in a drastic narrowing of the bundle of Sub(q → 0) asymptotics. Actually, the dependences Sub(q) have the same first q-derivative at q → 0. This fact has an important consequence and will be considered below in detail.
In our previous publication [59], we found that, in the case of a two-layer dispersionless systems with a given ratio ε 2 /ε 1 , all curves such as those presented in Figure 4a can be reduced to a single curve ε eff (z). It is remarkable that this unique curve could be used not only to calculate the r-dependences of the interaction between two charges with z = z , but also the z-dependence for a charge pair with z = z for a fixed r. Now, the availability of the second interface and the new length parameter L makes the existence of a unique function ε eff (z) impossible.
Thus, the calculation of the integral in Equation (3) by expanding the denominator in its integrand in the parameter θ makes it possible to calculate the interaction between two charges with an accuracy set in advance (in most cases, we required a relative accuracy of 10 −6 in series calculations). However, as was indicated above, the convergence rate of the series can be rather slow. The reason can be understood from Equation (11) and formulas for {F i } in Appendix C. In the case z = z , the functions f i (z, z ) are either comparable with 1 (if they depend on the difference z − z) or with the sum z + z = 2z. If r and z are small, both the factor (−θ) j and the number 2j in the denominator of Equation (11) have a substantial effect on the magnitude of the next series term. Otherwise, the effect of increasing 2j is insignificant, and the convergence rate is almost solely determined by the parameter θ. Thus, we can describe the domain of slow convergence of the series in Equation (11) as (r 1 and/or |z| 1) and |θ| → 1.
Hence, there arises a problem of constructing rather compact and accurate approximations that could be expressed in a finite analytical form.  Figure 6. The same as in Figure 5, but for symmetric structure E.

Effective-Exponential Approximation
The main difficulty in calculating the integral in Equation (3) is the existence of the denominator in Equation (5) in its integrand, which incorporates multiple image reflections [15,27,106] and comprises the interface polarization contribution to the electrostatic interaction energy. While obtaining an expression that would represent integral (3) rather well in the whole interval r ∈ [0, ∞), the most trivial, but reasonable, idea consists in that the denominator should be replaced by a function that can fit it well, but allows analytical integration. The class of such functions is not wide. The following reasons have to be taken into account.

•
The function Den −1 (q, ε) changes monotonically from (1 + θ) −1 at q = 0 to 1 at q = ∞, and, as discussed above, it contains no singularities. • Since Equation (3) involves the Bessel function in the integrand, it is very important for the new function to satisfy the required asymptotics both at q → 0 and q → ∞.
Therefore, we approximate the function Den −1 (q, ε) by the function Note that this approximation, similar to Den −1 (q, ε) itself, does not include the z-coordinates of the charges, so that it remains valid in the general case z = z . The parameter a is chosen to satisfy the equality The integrals on both sides are calculated analytically, and we arrive at the solution Substituting the function in Equation (23) into Equation (3), we obtain where a is determined by Equation (24). Below, this approximation will be referred to as E0. Note that the sum in Equation (25) has a finite number of summands. It is remarkable that the final result can be obtained in the analytical form for any combination of problem parameters, i.e., the dielectric constants of the layers, the width of the inner layer L, the charge locations Z and Z , and the lateral distance between the charges R.
As an example, let us derive the approximation E0 for the interaction between two charges with z = z located in the interlayer. For this purpose, we should use the set in Equation (20) for w 22 . Accordingly, Equation (25) reads where a is described by Equation (24), and θ by Equation (6). The other variants are calculated similarly. Hence, the proposed approximation of the charge-charge interaction in a three-layer structure reduces to a simple algebra. In Figure 7 (left), the results calculated using the series expansion (see Section 4) of Equation (3) and the analytical approximation in Equation (25) are shown for z = z and our toy structure B (the results for structures A and C are qualitatively similar). One can see that, at comparable ε i -values (corresponding to relatively moderate θ-magnitudes, see Table 1), our interpolation formula (the dashed curve) describes the behavior of the electrostatic interaction (the dotted curve) excellently within the whole r-interval [0, ∞) for a wide range of the parameter z. The relevant curves are almost indistinguishable from each other, so that the right panels demonstrate the relative deviation of the E0 approximation from the exact rw(r)dependence (the dashed curves). However, what about the domain of slow convergence in Equation (22)? Figures 8 and 9 illustrate the results of corresponding calculations for symmetric toy structures D and E, respectively, in which the dielectric permittivity ε 2 of the interlayer differs rather strongly from that of the covers (see Table 1). One can see that the approximation (the dashed curves) is not so good as in the previous case. To elucidate this issue, we should analyze the approximation Inv(q, ε) in more detail. Actually, Inv(q, ε) is a function of q and θ. This is also true for its parent Den −1 (q, ε). Figure 10a compares those two functions for several θs, and Figure 10b demonstrates their ratio in the whole θ-interval. As long as θ −0.5, the approximation is satisfactory. However, at θ −0.5, the function Inv(q, θ) may strongly deviate from Den −1 (q, θ), especially at small qs. In Appendix D, we show that, in the integral in Equation (3) with a large enough r, the Bessel function plays a role of a generalized δ-function, so that the interval of small qs makes the dominant contribution to this integral. In our opinion, just this situation emerges at θ −0.5, so that Inv(q, θ) contributions that are "overestimated" at small qs distort the true value of the integral. Moreover, the situation can be even more complicated. The summands in the sum in Equation (7)-actually, they are exponential functions-may strongly compensate one another at small qs.
In this case, a q-interval that is more distant from q = 0, where the function Inv(q, θ) is "underestimated" (see Figure 10b), can insert the dominant contribution. Anyway, the final result will be unsatisfactory. A receipt to overcome this difficulty was proposed in Section 3. In particular, we should integrate the leading term(s) analytically and regard the remainder as a correction (see Equation (13)). The latter should be approximated using, e.g., the method described in this section. More specifically, we used the procedure in Equation (13) to every ith term in the sum in Equation (7) for which f i (z, z) = 0. The extracted terms were combined into a single pure Coulombic component and extracted into the Coulomb tail, whereas the remainders, together with other F i terms, were treated in the framework of the effective-exponential approximation considered in this section. In essence, this is a version of the power grouping procedure described in Section 4. The dash-dotted curves in Figure 8 demonstrate the results obtained for symmetric toy structure F following this method. It is evident that the application of the latter leads to a substantial improvement of final results.
We emphasize that it is worth doing the term grouping procedure in any case. For instance, in the case when the charges are in medium 3 (see the set {F i } in Equation (21)), one can see that F 1 + F 4 = 1 ε 3 Den(q, θ), so that those terms produce a single and exact common contribution. At the same time, the terms F 2 and F 3 can be grouped. As a result, we obtain a single integral for the approximation.
If we consider the transformation of Equation (12) into Equation (13) more thoroughly, we understand that the grouping procedure applied to the terms with f i (z, z) = 0 is nothing but an extraction of the exact short-range (r → 0) asymptotics into the Coulomb tail. This approximation will be denoted as ES. However, in Figure 8, we had troubles with large r-distances. Therefore, there arises an idea to preliminarily extract the exact long-range (r → ∞) asymptotics. It is easy to understand that the corresponding coefficient equals Sub(q = 0). This approximation is denoted as EL, and its results are shown in Figure 8 by solid curves. Note that the ES and EL approximations coincide in the symmetric case if |z| > 1 2 , because the short-and long-range asymptotics are identical for the electrodes. The difference between those approximations in the whole interval of z-variation can be observed in Figure 7.
The same extraction procedure of different Coulombic asymptotics was also applied to symmetric toy structure G. The results obtained were found to be similarly satisfactory (see Figure 9). Thus, the extraction of either short-or long-range asymptotics improves, as a rule, the approximation accuracy. Nevertheless, the comparison of Figures 8 and 9 testifies that the selection between those asymptotics must be done carefully. Again, no recommendation can be made a priori.
To summarize this section, the effective-exponential approximation is very simple and seems to be very useful and easily applicable in practice for estimating the behavior and magnitude of the interaction energy w(z = z , r) in heterostructures.

Long-Range Asymptotic Power Sums
The easiest way to fill the lacuna in the effective-exponential method is to construct a long-range asymptotics for the integral in Equation (3) in the form of asymptotic sums. Really, since the function Sub(q, z, z ; ε) does not depend on r, this integral belongs to the class of integrals in Equation (A33) considered in Appendix D. Therefore, we may evaluate its long-range (r → ∞) asymptotics using the properties of J 0 (qr) as a generalized δ-function. This means that only the properties of the function Sub(q, z, z ; ε) at small qs are relevant. Hence, we should first expand the Subintegrand into a Taylor power series in q, Then, every term in the expansion should be integrated with J 0 (qr) using equations of Appendix D. As a result, we obtain an asymptotic series of the type which contains only odd powers of r. The coefficients as in the expansion in Equation (28) depend on the problem parameter set, which, in our case, includes the dielectric permittivities of the media and the (z, z )-coordinates of the charges, and does not depend (at i ≤ N a ) on N a . The analytical results obtained for a i can be obtained straightforwardly, but they are rather cumbersome. Therefore, as an example, we present the approximation formula for w(z = z , r; ε) only in the case when the both charges are located in the interlayer (− 1 2 ≤ z = z ≤ 1 2 ), and retain only the first two terms: i = 1 and 3 (this approximation are denoted as T3): where a 1 = 2 and The complexity of expressions for the higher-order coefficients increases drastically with their order numbers. Nevertheless, they remain to be ordinary algebraic expressions similar to those presented above.
One should note that the coefficient a 1 is z-independent, as it has to be (see Equation (16)), whereas the coefficient a 3 is a quadratic polynomial of z including the non-zero linear z-term. However, in the symmetric geometry (ε 1 and the linear z-term disappears, which is important for further consideration. In Section 5, we emphasize an unambiguous relation between the Sub(q) and w(r) functions. In particular, this concerns a relation between the short-wave asymptotics Sub(q → 0) and the long-range asymptotics w(r → ∞) (see the relevant panels in Figure 4; in the case of symmetric structures, see Figures 5  and 6). Here is a corresponding short-wave counterpart of the long-rang asymptotics in Equation (29): In the symmetric geometry (ε 1 A correlation between those two sets of formulas is obvious. The latter set explains the collapse of the divergent bunch of the Sub(q) curves for various zs in the case of non-symmetric structures (see Figure 4c) into a narrow bunch of the corresponding curves divergent only in the second derivative in the case of symmetric structures (see Figures 5c and 6c).
In Figure 11, the exact long-range w(r) profiles (dotted curves) calculated for toy structure A are compared with their asymptotic sums (Equation (28)) calculated for N a = 3, 5, 7 at certain characteristic z-values. We selected this structure with comparable εs to demonstrate that the long-range asymptotics of w really depend on z, although weakly. One can see that the lower end of the convergence interval for the sums in Equation (28) becomes larger as the charges go farther from the slab (the parameter z increases by magnitude). It is also evident that the further (N a > 3) expansion of the series in Equation (29) is not required for qualitative calculations. However, again, what about the domain of slow convergence of the regular method, where the effective-exponential approximation is not efficient? Figure 12 demonstrates that the method of long-range asymptotic power sums resolves, to some extent, this problem. However, the corresponding interval where the sum approximates the dependence w(r) can be not enough to cover the interval where the effective-exponential method works badly. Figures 11 and 12 also illustrate the following fact. For any finite N a , the sum in Equation (28) is divergent at r → 0. Therefore, the w(r) profile is very sensitive to the values of the expansion coefficients a i . In the case of comparable ε i s, the change of z modifies those coefficients quite enough for this modification to appreciably affect the sum profiles in Figure 11. This fact confirms the validity of Equation (32). In the case ε ε , the z-dependent term prevails in the multiplier a 3 , and small variations in the long-range behavior of the dependence w(z) can be observed (Figure 12a). However, if ε ε , the validity of the expansion in Equation (29) is universal but less satisfactory than in the opposite case ε ε (Figure 12b). Attention should be paid that the asymptotic series in Equation (29), as well as the omitted terms in it, was obtained directly from the integral in Equation (3) without any approximations. Therefore, this series is exact and can serve as a reference to verify the quality of various approximate expressions, at least in the long-range limit.

Preliminary Remarks
The Rytova-Keldysh approximation (RKA) for the electrostatic interaction in three-layered structures was for the first time suggested and analyzed by Rytova while studying the screened potential of a point charge in a slab [113] and later on by Keldysh (in a slightly different approach) for the charge-charge interaction in the same geometry [48]. However, both authors analyzed the interaction between two charges located at different sides of the slab, so that the z-coordinates of the charges were fixed (z = − 1 2 and z = 1 2 in our notation), and the interaction energy w was z-independent. The concrete results obtained in those articles [48,113] are often applied to a number of charge configurations without specifying them and applying the RKA formula mostly to exciton physics in heterostructures [68,72,73,[76][77][78]98,101,[114][115][116][117][118][119].
In Section 9, we consider Coulomb interaction in the configuration z = − 1 2 and z = 1 2 in more detail. At the same time, in this section, we analyze how the RKA can be modified to include the z-coordinates of the charges into consideration, which is necessary, in particular, to improve the numerical results. In this connection, we would like to emphasize that under the term "Rytova-Keldysh approximation", we mean an analytical approximate method used to calculate integrals such as that in Equation (3)  The following circumstances have to be taken into account. In essence, the RKA was developed to calculate the long-range asymptotics of the integral in Equation (3). It is based on the property of the Bessel function J 0 (x) to play the role of a generalized δ-function (see Appendix D). Then, if we are interested in the long-range asymptotics of the integral in Equation (3), the task of approximating the profiles of Subintegrand Sub(q) in the whole interval 0 ≤ q < ∞, which are rather complicated (see, e.g., the relevant panels in Figures 4-6), becomes not urgent. Now, we have to choose such an approximation for Sub(q) that (i) the resulting formula be integrable analytically, and (ii) the approximation be accurate enough at small qs.
Formally, the RKA is a two-parametric approximation of the function Sub(q) in the vicinity of q = 0 selected in the form Hence, we need two independent quantities (conditions) to determine the parameters A and D. Among such quantities, the most natural are the values of the function Sub(q = 0) and its derivative d dq Sub(q = 0). All other possible sets of two independent parameters must be related to this one. However, for a system with a given set {ε i }, the value (see Equation (33)) Sub(q = 0) = 2 ε 1 + ε 3 (37) and, thus, the ratio A/D, is fixed (see Equation (15)). Therefore, the derivative value d dq Sub(q = 0) acquires an extremely important significance.
The approximation in Equation (36) is applicable only if D > 0. In this case, Equation (3) is integrable analytically (see Equation (A41)), and the approximation formula for w(r) looks like (see Equation (A41)) It has a logarithmic singularity at r → 0 instead of the proper pole one [68,101]. Therefore, it can never be applied to calculate the short-range w(r)-asymptotics, although the r-interval of its applicability can be rather wide for some {ε i } sets and z-values.

RK0 Approximation and Its Modifications
In the basic RKA version (we refer to it as RK0), the numerator Num(q) in Sub(q) (see Equation (4)) is substituted by its value at q = 0, whereas the denominator Den(q) is approximated by a linearized Taylor series in q. After the corresponding normalization of the numerator and denominator, we arrive at Equation (36), where These parameters ensure both the correct value of the fraction in Equation (36) at q → 0 and its vanishing at q → ∞. One should note that all information about the position of charges (z and z ) was contained in the numerator Num(q), and hence it is lost altogether in the RK0 approximation. As was indicated above, this assumption was sufficient for the original approximation to be applicable in the across-slab configuration of charges in the layered system [48,113] (see also Section 9).
For the right hand side in Equation (36) to have no singularity at q ≥ 0 (i.e., D > 0), the parameter θ must be negative, i.e., −1 < θ < 0. This condition is realized if ε 2 > max(ε 1 , ε 3 ) or ε 2 < min(ε 1 , ε 3 ), so that toy structure A becomes excluded from consideration at once. However, this condition is always valid for symmetric (ε 1 = ε 3 ) heterostructures. In any case, since the RK0 approximation is a monotonously decreasing q-function (therefore, its derivative at q = 0 is always negative), it will result in a monotonously increasing r-asymptotics.
For the other toy structure C, the corresponding dashed curve in Figure 13a demonstrates that the RK0 approximation is too rough to provide reliable results except for very large r, which is no wonder, because RK0 was constructed to describe this limit. As is evident from Equation (39), both parameters A and D are fixed, so that the RK0 profile is also constant and cannot trace the variation of the rw(r)-profile with z (the dotted curve). The situation is even worse for toy structure B (Figure 13b), because the RK0 and rw profiles disperse in different directions from the asymptotic point. Indeed, the actual derivative Hence, we have to modify Equation (36) and approximate the derivative d dq Sub(q = 0) in such a way that it would improve the fitting of the original Sub(q) dependence. In so doing, we are not allowed to change the linear binomial in the denominator, because more complicated functions will make the integral in Equation (3) non-integrable analytically. Instead, we may try to expand the numerator Num(q) (see Equation (4)) into a Taylor series. In the case of linear expansion, it is easy to see that, after simplification, the approximating function can be now written in the form (the approximation RK1) where D is the same as in the approximation in Equation (36), i.e., it is given by Equation (39), but A and B are determined on the basis of the both values Sub(q = 0) and d dq Sub(q = 0). The variation of the parameter B corresponds to the shift of dependence (36) along the ordinate axis, so that now we can reproduce dependences with positive d dq Sub(q = 0): for this purpose, A must be negative and B positive and large enough. The actual values of the parameters A and B in the geometry z = z are as follows: One can see that, in the symmetric configuration (ε 1 = ε 3 ), the both parameters become z-independent. Of course, the latter fact follows from the z-independence of the derivative d dq Sub(q = 0). This is illustrated by Figures 5c and 6c, as well as by Equation (32). The latter also testifies that we have to expand the numerator Num(q) up to the second-order term to ensure that the approximate Sub(q) and, hence, rw(r) depend on z. After a straightforward simplification, we obtain (the approximation RK2) The parameter D, as expected, remains the same, but new expressions for the parameters A, B, and C are cumbersome, and their explicit formulas are omitted here. Every term in the RK2 approximation multiplied by the Bessel function J 0 (qr) is integrable.
The following circumstance is of interest. The coefficients in the series expansion of Num(q) for given z and z do not depend on the expansion order. However, when Sub(q) is transformed to the forms containing the RK term (Equation (36)), the parameters, firstly A and then B, change from one approximation to the other. The expressions for A and B in the approximation in Equation (42) differ strongly from Equation (41) for the approximation in Equation (40). At the same time (see Equation (A37)), the last term in Equation (42) does not contribute by no means to the rw-dependence, although the Subintegrand approximation does change! It looks as if the role of this term consists in a proper renormalization of the coefficients in the approximation in Equation (40).
All the aforesaid is illustrated in Figure 13. One can see that the account of the next terms in the series expansion of Num(q) strongly improves the approximation quality. This is especially true for toy structure B (Figure 13b), for which the behavior of the basic RK0 approximation is totally wrong.

Original RKA
Let us pay attention that Equation (36) has an extremely felicitous property. If the {ε i }-driven parameter θ tends to −1, the singularity point q = θ+1 2θ < 0 approaches zero. As a result, the contribution to the integral in Equation (3) from the vicinity of the left-end point q = 0 increases and, at θs sufficiently close to −1, becomes dominating. The limit θ → −1 is realized in two cases: (i) ε 2 min(ε 1 , ε 3 ); and (ii) ε 2 max(ε 1 , ε 3 ). As noted in the previous section (see also Equation (34)), the derivative d dq Sub(q = 0) is negative in the former case, so that it must be excluded from consideration. Let us take symmetric structures (ε 1 = ε 3 = ε , ε 2 = ε ) to illustrate this scenario. In this case, the parameters A and D in Equation (36) look like We recall that the results of the approximation in Equation (36), which was used to calculate long-range asymptotics, remain to be z-independent. This model is satisfactory at large ε /ε ratios. Indeed, even for toy structure D with {ε i } = 1 : 1.5 : 1, the profile of the dependence Sub(q) near q = 0 changes weakly for various zs [see Figure 5c; note that the curves in that figure span the z-coordinate values |z| ≤ 10, whereas the slab is located at |z| ≤ 1 2 ]. In the case of toy structure F with {ε i } = 1 : 10 : 1, this change is much weaker according to Equation (35). Therefore, the calculations in this section can be considered as those made for the charges located, e.g., near the middle of the slab (z = z ≈ 0). Figure 14 illustrates the quality of the RK0 approximation for various values of the ratio ε /ε > 1 (cf. the dotted and dashed curves). A drastic improvement of the RK0 approximation with the increase of ε /ε is obvious. It is worth noting that this improvement takes place in a wide interval of r-values rather than only at r 1. Nevertheless, as was indicated above and is illustrated in the inset of the bottom panel in Figure 14, the difference between the logarithmic (appropriate to RKA) and pole (genuine Coulombic) singularities at r → 0 survives at any ε /ε -ratio, so that this approximation cannot be applied if one is interested in the electrostatic interaction behavior at very short distances between the charges. Thus, one can see that two factors-large θ-magnitudes (θ → −1) and the property of the Bessel function J 0 (qr) to play the role of δ(q)-function in integrals of the type in Equation (3) at large rs-"work cooperatively" to achieve the same goal: to allow the close vicinity of q = 0, where substitution of Equation (36) is justified, to make a dominant contribution to the integral in Equation (3). If either of those conditions is not met, the other can compensate the absence of its counterpart by becoming stronger.
Keldysh [48] made a next step and, in accordance with the relationship ε 2 max(ε 1 , ε 3 ), coarsened the parameters A and D (Equation (39)): or, in the symmetric case, We refer to this approximation as RKorg. Amazingly, but in most cases, the "coarsened" RKorg approximation turns out to be more successful than the seemingly "more accurate" RK0 one (see the solid curves in Figure 14)! Moreover, this is also true for non-symmetric structures with comparable ε i s. In a good many cases, the RKorg approximation turns out to be even more accurate than the corresponding RK2 one (see Figure 13a). Of course, this is true only if the RKorg approximation is applicable (see Figure 13b).
To explain this fact at the present stage of discussion, let us recall that the quality of the approximation in Equation (36) is determined not least by its derivative value at q = 0. In Figure 15, the ratios of the derivatives d dq RK0(q = 0) and d dq RKorg(q = 0) to the actual derivative d dq Sub(q = 0) are compared with each other. There is no doubt that the RKorg approximation wins the competition.

Moreover, let us rewrite the RKorg approximation in the form
whence it becomes evident that it corresponds to a fractionally rational (homographic) approximation to Num(q). The capabilities of functions belonging to this class strongly exceed the capabilities of Taylor series. In particular, fractionally rational functions are used to approximate complicated special functions to a high accuracy within rather wide intervals (see, e.g., [120,121]).

Fitting Approximations
The success of the RKorg approximation with the coarsened parameters in Equation (44) or Equation (45) brought us to the following idea. In essence, the coarsening procedure detached the original physical meaning from the quantities A and D. We propose to go further. Namely, let us declare them fitting parameters and use the values of both Sub(q = 0) and d dq Sub(q = 0) to determine them. As an example, we obtain (the approximation F0) for charges located in the slab (− 1 2 < z < 1 2 ), and if they are in medium 3 (z ≥ 1 2 ). The analytical formulas for the coefficients in the models in Equations (40) (the approximation F1) and (42) (the approximation F2) are very cumbersome algebraic expressions and are not shown. Attention should be paid that now every approximation requires an extra condition, i.e., additional information about Sub(q). In particular, for the approximation F2, it is necessary to know the expressions for Sub(q = 0) and d i dq i Sub(q = 0) with i = 1 to 3. One can see that the fitted values depend on z even in the lowest approximation, but (as expected) only for non-symmetric structures. Since the denominator of the parameter D depends on z, it, in principle, can change its sign. Therefore, the sought approximation may exist only in definite z-intervals. At the same time, the emerged z-dependence of the coefficients A and D can make the fitting F0, F1, and F2 approximations applicable for those {ε i } sets (again, within definite z-intervals), where their non-fitting counterparts fail. Figure 16 demonstrates the results of calculations for toy structures A, B, and C with comparable ε i s. The results testify to a good quality of fitting approximations, which allows them to be used down to very short distances between the charges. The results obtained for toy structures D and E are even better and, hence, are not displayed. At the same time, there are no F0 plots in some panels, which means that this model failed for the corresponding parameter sets. Now, we can explain the success of the RKorg approximation from a more general viewpoint. This approximation is an intermediate point on the way from the initial, primitive RK0 approximation to its more complicated fitting extension. Keldysh made a partial correction, which is, generally speaking, necessary for the influence of the Num(q) dependence on the overall behavior of w(r) to be taken into account. Maybe he was lucky to guess "coarsening" approximation (44) leading to the description improvement, because the coarsening corrects the "theoretically accurate" coefficients in Equation (39) in the proper direction.

Remarks on the Margin
We have paid so much attention to the RKA, because it is used rather often in recent studies of heterostructures [73,98,101,[114][115][116][117][118][119], although sometimes beyond the scope of its applicability [122]. This approximation was considered as an auxiliary tool with respect to the effective-exponential approximation (Section 6) in the domains of problem parameters, where the latter approximation gives unsatisfactory results. Both approaches can be used in combination to check the results obtained in the framework of other approximations. Their compatibility follows from the complete equivalence of their long-range asymptotics (irrespective of the RKA versions, except for the RKorg with its coarsened parameters (Equation (45))) obtained to an arbitrary order, both with each other and with those obtained in the exact method of long-range asymptotic power sums (Section 7). However, the RKA is rather cumbersome, and totally fails in the very important case of heterostructures with close values of dielectric constants [62,63,66,[108][109][110]122]. Therefore, the effective-exponential approximation is preferable, as a rule. Moreover, one should never use any version of the Rytova-Keldysh approximation to study the interaction of charges at short distances, r ≤ 1.

Preliminary Remarks and Analytical Formulas
For many fundamental and practical problems [65,66,73,84,123], it is important to elucidate how two point-like charges interact with each other if they are located at different interfaces of the slab in the three-layer system. The interaction of charges through a slab is of special interest for the problems of exciton spectra in heterostructures [48,[67][68][69][70][71][72][73][74][75]77,78,119] as well as those concerning exciton (electron-hole) superfluidity in layered systems [46,47,[79][80][81][82][83][84][85][86] and related topics [75,[123][124][125][126]. In some of those works, the actual electron-hole interaction was considered to be a conventional Coulombic attraction, whatever the relationship between the dielectric permittivities of the media involved. Meanwhile, any proper analysis of the electron-hole pairing must include a detailed account of the electrostatic interaction screening by all three layers. That is why a number of attempts were made to describe, more or less correctly, the charge-charge interaction as applied to the indicated phenomena [11,48,70,113]. However, in the problems concerned, the RKorg approximation (Equation (45)) was and remains the most popular theoretical foundation of calculations.
Broadly speaking, the spatial dispersion of ε i (k) should be taken into account, especially if itinerant charge carriers are present. Of course, the corresponding results will be cumbersome and less transparent. Therefore, since the problem has not been analyzed even in the adopted here simplest possible CDP model, we believe that it is worthwhile obtaining an insight in the framework of the classic electrostatic approach. Hence, a firm basis should be developed in order to check the validity of the existing simplifications in the framework of the theoretical scheme (Equations (A1)-(A24)).
As indicated above, the required {F i } set can be obtained from any of the sets for w 12 , w 13 , w 22 , and w 23 (see Appendix C) by putting z = − 1 2 and z = 1 2 in them. Let us use the set in Equation (A29) as the simplest of them (it includes only one term). Substituting Equation (A29) into Equation (3), we see that the interaction of charges located in the both outer layers is "the true Coulombic one" only if θ = 0, i.e., if ε 2 equals either ε 1 or ε 3 , or, in other words, if the system is two-layer (or one-layer if ε 1 = ε 2 = ε 3 ). When θ = 0, we have Below, we consider the most interesting, from the practical viewpoint, symmetric case (ε 1 = ε 3 ≡ ε = ε 2 = ε ). Then, the parameter θ becomes and, together with the denominator in Equation (5), does not depend on the interchange ε ↔ ε . As a result (see Equation (3)), we obtain an interesting symmetry relation for the interaction between two point charges located in different cover layers of the studied heterostructure: Of course, this relation is valid in the limiting case z = − 1 2 , z = 1 2 as well. In particular, Equation (3) acquires the form which is evidently consistent with the relationship in Equation (51). The integral I in Equation (52) can be calculated analytically if the charges are arranged vertically (r = 0). The result is as follows: In the symmetric case, which is considered in this section, θ < 0, and we obtain

Results of Calculations
Let us proceed from the calculation of Subintegrand. According to Equation (49), it is equal to Its Taylor series up to the linear term equals One can see that the derivative d dq Sub(q = 0) is always negative, so that the RK approximation is also valid, at least formally. Note that Keldysh overlooked this possibility, having assumed that his approximation is valid only if ε 2 > max(ε 1 , ε 3 ). We emphasize that the last statement concerns only the specific configuration of charges located at different interfaces.
In Figure 17, the dependences Sub(q) are shown for all toy structures. The analysis of the curves presented in this figure and the calculation results for the w(r) dependences that are discussed below completely confirm conclusions made in the previous sections on the validity domain of each examined approximation. In particular, this concerns the mutually complementary character of the effective-exponential and RK approximations. From the viewpoint of their applicability, the former works excellently in a wide range of problem parameters, except for large r and incomparable ε values in tandem; the latter works as if it were specially designed to fill this lacuna.
The curves in Figure 17 look very similar: a monotonic, quite rapid (exponential) reduction to zero at q → ∞. Therefore, we may confine the analysis only to the symmetric configuration. Moreover, since in this case, the relationship in Equation (51) takes place, we consider only configurations with ε 2 > ε 1 = ε 3 .  Figure 18 illustrates the results of calculations for toy structure D ({ε i } = 1 : 1.5 : 1) using the regular series expansion method expounded in Section 3 (the dotted curve) and the E0 approximation in Equation (55) (solid curve). One can hardly detect the difference between the curves, because the relative approximation error does not exceed 1.4 × 10 −4 in the whole interval 0 ≤ r < ∞. The extraction of the long-range asymptotics (the approximation EL) allows this error to be reduced to 5.9 × 10 −5 . Note that the isolation of the short-range asymptotics (the approximation ES) has no sense in this case, because the charges do not approach each other closer than the slab width.  Figure 18. Dependence sw(r) and its effective-exponential approximation E0 for symmetric structure D.
Owing to the relationship in Equation (51), the w(r) profile for toy structure E ({ε i } = 1.5 : 1 : 1.5) can be obtained from Figure 18 by changing the ordinate axis scale. The approximation error values remain the same.
The picture becomes more complicated if the difference between the dielectric constants in the slab and covers is large. As indicated in Section 6, the effective-exponential approximation cannot provide sufficient accuracy at large inter-charge distances in this case and should be appended, e.g., by the RKorg approximation. In Figure 19a, we demonstrate such a combination for toy structure F ({ε i } = 1 : 10 : 1). An even more pronounced example of this "cooperation" is illustrated in Figure 19b for toy structure H with a very large difference between the dielectric constants in the slab and covers ({ε i } = 1 : 100 : 1). One can see that the seemingly "more accurate" EL approximation can produce a physically unreasonable result (the negative static dielectric constant). At the same time, the E0 and RKorg approximations provide an almost smooth approximation of the rw(r) dependence with a high accuracy in the whole approximation interval 0 ≤ r < ∞.
Our calculations clearly demonstrate that, e.g., in the exciton physics of heterostructures [48,[67][68][69][70][71][72][73][74][75][76][77][78]98,101,[114][115][116][117]119], where the dielectric permittivities of the media are usually of comparable values, it is better to use our effective-exponential approach rather than any version of the RKA. Nevertheless, the latter might turn out successful for some parameter sets, although those approximations are not at all universal (in particular, they contain an inevitable logarithmic singularity at small arguments). Furthermore, the effective-exponential approximation is attractive owing to its simplicity and convenience, because the specific approximate expressions for w(s) are algebraic ones.

Conclusions
In this work, we calculated the charge-charge electrostatic interaction W in a three-layer structure with plane interfaces. Exact formulas for the dependence of W on the relative lateral, R, and perpendicular with respect to the interfaces, Z and Z , coordinates of the charges were presented. In the framework of the simplest approach, where the dielectric functions of the media were treated as dispersionless constants, simple approximations allowing finite analytical expressions were derived and discussed. They allow one to obtain an insight into the problem, as well as to use them in practical applications, at least for estimation purposes. The results obtained can serve as reference ones when developing the theory further. The results found in the framework of different approximations were compared and analyzed.
Two specific charge arrangements were analyzed in detail. This is the case when the Z-coordinates of the charges are identical (Z = Z ) and when the charges are located at different interfaces. It was shown that the Rytova-Keldysh approximation (RKA) [48,113], which is widely used, e.g., in the physics of excitons [73,74,76,101,114,115,125,129,130], turns out to be unsatisfactory in heterostructures with similar dielectric constants of the media involved. Instead, we proposed a new analytical approximate formula (the effective-exponential approximation) that fits well the exact integrals [6,11,15] in a wide range of their parameters. The spectrum of applications of the proposed theory can be rather wide. For instance, it can be used in the studies of possible electron-hole superfluidity in layered structures [46,47,[79][80][81][82][83][84][85][86]. media 1 (Z < 0) and 3 (Z > L) separated by a homogeneous media 2 (0 < Z < L)-each of the ith medium (i = 1, 2, 3) is characterized by the dielectric permittivity function ε i (k, ω)-can be found using several equivalent approaches [6,11,14,15] based on the assumptions of the infinite interface barriers for particles forming the media involved, the specular reflection of the charge carriers at the corresponding interfaces, and the absence of relativistic retardation effects that take into account the finiteness of the light velocity [131][132][133][134][135]. In particular, the energy W QQ (S) can be written in the form where the subscripts i and j denote the media where the charges Q and Q , respectively, are located; R is the horizontal distance between the charges; Z and Z are their vertical coordinates; and J 0 (x) is the Bessel function of the first kind. The functions where δ ij is the Kronecker delta, are called Green's functions of the three-layer Coulomb problem [1]. They describe the screened electrostatic interaction and allow both the charge-charge interaction and the self-energy of the charge in the layered system (the image charge energy [6,11,15]) to be calculated. In the reference frame described above, the auxiliary quantities entering Equation (A2) are as follows: where

Approximation of Dispersionless Dielectric Permittivities
Equations (A1)-(A21) allow one to obtain electrostatic interactions between arbitrarily located charges in a three-layer structure made of arbitrary materials. Hereafter, in accordance with the aforesaid, we confine the consideration to the simplest possible case of the dispersionless ε i (q, k ⊥ , ω) ≡ ε i = const in all three media. This approximation is valid not only for the universal "medium" -vacuum-but also for insulators and semiconductors with large band gaps. Although the approximation concerned seems trivial at first glance, we noticed that even this case has not been examined in the literature in full. In this approximation, we arrive at the three-constant (ε i ) parameterization of the problem, and the quantities a 1,S,A,3 (Z) can be integrated analytically: If the argument Z goes beyond the indicated limits, the functions a S,A (Z) are recalculated according to the formulas a S (Z) = a S (Z ± L), a A (Z) = −a A (Z ± L). (A26)

Appendix B. Parameter θ
The parameter plays a key role in the described theory. In particular, it governs the convergence rate of the series used to calculate the dependence w(r). Besides that, it is responsible for the applicability of the RK approximation for w(r) at r → ∞.
Since in the framework of the classical electrostatics [106] all the dielectric constants entering the block in Equation (6) are positive (we note that it is not at all the case when the spatial dispersion of the dielectric permittivity is taken into account [136]), it is convenient to evaluate θ in terms of the relative differences between the relevant arguments. For instance, for the quantities x and y, the relative difference equals δ(x, y) = x − y x + y .
It varies from −1 at x y to +1 at x y and equals zero at x = y. In some sense, this function is inverse to the Kronecker delta δ ij : the latter equals unity if its arguments are equal to each other, and becomes zero otherwise. Therefore, we use a bar in the notation to emphasize this difference.

Appendix C. Constituents of Interaction Energy w(z, z )
In this appendix, we present full sets of summands in the sum in Equation (7) that are required to calculate the interaction energy between two point-like charges in the three-layer system concerned. In particular, these are (we use the notationδ(x, y) for the relative difference, which was defined in Appendix B): • for w 11 (z ≤ z ≤ − 1 2 ): N = 4, • for w 12 (z ≤ − 1 2 ≤ z ≤ 1 2 ): N = 2, • for w 13 (z ≤ − 1 2 , z ≥ 1 2 ): N = 1, (ε 2 , ε 1 ) exp −q z + z + 1 , • for w 23 (−1/2 ≤ z ≤ 1/2 ≤ z ): N = 2, (ε 2 , ε 1 ) exp −q z + z + 1 , To avoid ambiguities but without any loss of generality, we confine ourselves to the case where the coordinate z of the charge in medium i does not exceed the coordinate z of the charge in medium j. It is evident that, due to the symmetry of the problem itself and the reference frame, the set for w 23 can be derived from the set for w 12 and the set for w 33 from the set for w 11 if one makes the obvious changes of indices.
All sets match in the corresponding limiting cases. For instance, the set for two charges that are located at different interfaces can be obtained from the sets for w 12 , w 13 , and w 22 by putting z → −1/2 and z → 1/2. Summation of the obtained terms brings us to the same result in each case.
where {a} stands for a certain parameter set (it does not include the coordinate r), and G is a smooth enough and finite function, so that this integral converges. However, the set of G-functions for which the integral in Equation (A33) can be calculated analytically is rather restricted [137]. The quality of the approximate calculation of this integral is associated with the behavior of the Bessel function J 0 (x), which is illustrated in Figure A1. The lower panel demonstrates that, if we take a relative accuracy of, e.g., ±10% and G(q, {a}) = const, the main contribution to the integral is given by the interval rq ≈ [0, 60]. At r = 1, this interval corresponds to the interval q = [0, 60] in the integral in Equation (A1). Extending the upper integration limit beyond 60 will only provide a more accurate value within the selected error interval. Now, if the parameter r increases, say to 10, the interval rq = [0, 60] will correspond to q = [0, 6] in the integral in Equation (A1). In this sense, the Bessel function J 0 (qr) at r → ∞ behaves as a generalized δ-function, making relevant only the dependence of the function G(q, {a}) at small qs. Of course, the estimations made above are rough. They strongly depend on the behavior of the function G(q, {A}). In particular, at G(q → ∞, {A}) → 0, the role of J 0 (qr) as a generalized δ-function will be more effective. At the same time, if G(q, {A}) vanishes at q → 0, the upper integration limit has to be selected as large enough. Anyway, if the integral in Equation (A33) converges, we may calculate its long-range (r → ∞) asymptotics taking into account the behavior of the function G(q, {a}) in the vicinity of q = 0.
In this work, the formula [138] I(r, n, a) = ∞ 0 dqJ 0 (qr)q n exp(−aq) = (−1) n d n da n 1 √ r 2 + a 2 (A34) is actively used. In particular, at n = 0, we have and, at a = 0, the functions for various ns are: Note, that all values of I(r, n, a) with odd ns equal zero. Another well-known integral used in this work looks like [138] ∞ 0 dq J 0 (qr) where the both real-valued parameters r and a are positive (r > 0 and a > 0), and H 0 and Y 0 are the Struve and Neumann functions of the zeroth order, respectively. It is remarkable that similar expressions are typical of the polarization energy in the interface region if the spatial dispersion of dielectric permittivities is accounted for [139,140].