Correspondence between Electrostatics and Contact Mechanics with Further Results in Equilibrium Charge Distributions

It is common in mesoscopic systems to find instances where several charges interact among themselves. These particles are usually confined by an external potential that shapes the symmetry of the equilibrium charge configuration. In the case of classical charges moving on a plane and repelling each other via the Coulomb potential, they possess a ground state à la Thomson or Wigner crystal. As the number of particles N increases, the number of local minima grows exponentially and direct or heuristic optimization methods become prohibitively costly. Therefore the only feasible approximation to the problem is to treat the system in the continuum limit. Since the underlying framework is provided by potential theory, we shall by‐pass the corresponding mathematical formalism and list the most common cases found in the literature. Then we prove a (albeit known) mathematical correspondence that will enable us to re‐discover analytical results in electrostatics. In doing so, we shall provide different methods for finding the equilibrium surface density of charges, analytical and numerical. Additionally, new systems of confined charges in three‐dimensional surfaces will be under scrutiny. Finally, we shall highlight exact results regarding a modified power‐law Coulomb potential in the d‐dimensional ball, thus generalizing the existing literature.


DOI: 10.1002/andp.202300269
Thomson's model, [1] witnessing a paradigm shift from classical to quantum physics.In modern terminology, determining the minimum energy configuration of N identical point charges interacting via the Coulomb force, constrained on a 2D curved surface (mostly the S 2 sphere), is known as the Thomson problem.Its simplicity and geometric appeal has drawn the attention of mathematicians trying to find the optimal configurations. [2]The quantum mechanics treatment of low-density electron systems by Wigner [3] merely 30 years later predicted the existence of electron crystals, which added more theoretical ground to be explored experimentally.Both models play a major role in our understanding of equilibrium configurations of interacting particles in the cases of soft or hard confinements.
[6][7] Systems where electrons interact with the long range Coulomb force have attracted much attention because of the novel physics that emerges in low dimensions.The Thomson problem has provided interesting insights into quantum dots and Bose-Einstein condensates, [8] topological defects, [9,10] and charge-stabilized emulsions. [11,12]ypical examples are a 2D system of electrons realized in nanoscale semiconductor quantum dots [8] and electrons trapped on the surface of liquid He, [13][14][15] as well as the corresponding simulations. [16,17]n particular, electrons above a liquid He surface form a 2D electron gas that has a very low density and as such can be considered as a classical system of point charges where the Coulomb interaction between the electrons is the dominant effect.It is expected that, under these conditions, the electrons will form a 2D Wigner crystal state. [3]Let us recall that the Wigner crystal is a result of the interaction between the electrons and ripplons at the surface of liquid He, and the Coulomb interaction among the electrons.Thus, the Wigner crystal already involves elasticity of the He film.Let us bear in mind that the mapping of this interaction on the elasticity-theory problem is therefore somewhat more complex than a mapping of a pure electrostatic problem on the problem of the elasticity theory.A proper treatment can be found in Klimin et al. [18] It is known there that the minimum energy state of a classical infinite system of point charges corresponds to a triangular lattice, and the bulk properties of such a phase are fairly well understood. [19]However, the same conclusion cannot be drawn when dealing with finite systems of electrons.Advances in nanoscience and nanotechnology have led to the fabrication of devices consisting of a very few electrons.These devices are typically based on a 2D degenerate electron gas that can be controlled by patterned electrostatic gates.[22][23] The ensuing geometry and nature of the confinement potential is of paramount importance when it comes to describing the equilibrium state of finite and bulk systems of electrons.26][27] Among those confining potentials V(r), the hard-wall one represents a popular choice to study finite trapped systems of electrons.For instance, phase transitions in a finite 2D classical system of charged particles which are confined by a hardwall (alternatively, circular parabolic) well have been studied by Monte Carlo simulation techniques in the case of circular symmetry. [19,28]The effect of geometry on the configuration of particles and the melting process of a system of classical charged particles trapped by a circular hard-wall or by a parabolic external potential have also been studied. [29]Other confinements lead to vortices in superconducting films. [30,31]onte Carlo simulation methods to find the minimum energy configuration of a finite system of N interacting charges in a given confinement geometry is often an impossible task as N → ∞.Thus, the potential calculations of such systems in the continuum limit is often the only tool at hand.Let us recall that the equilibrium charge density for closed surfaces is determined by the curvature of the surface, and in the case of open surfaces (polygons, for instance), the edges start to play an additional role.As the number of particles grows, it becomes difficult to find the minimum energy of the system.Let us consider as an example the common case of charges bound to move in a circular confinement.The energy of a cluster of N charges confined to a disk of radius R, say by a hard-wall potential is given by (from now on in units of k e ≡ 1 4 = e = 1, unless otherwise stated) where As mentioned previously, finding the global minimum for a function such as E H is a very difficult task.The number of metastable states proliferate exponentially with N; consequentially the global minimum is obscured by a vast number of local minima with energies close to that of the global minimum.There exists a number of heuristic methods for such problems, such as quenched molecular dynamics or simulated annealing.Although there is no guarantee of finding the global minimum, it is possible to find states close to it.Systems involving a few thousand particles are computationally very much demanding, thus it is imperative to find alternative strategies.Therefore the first approach to the equilibrium of a system of N interacting charges via the usual Coulomb potential is provided by the continuum case, even if the system is of quantum origin.The standard theory that deals with the aforementioned confinements is given by potential theory.However, not many researchers are acquainted with such formalism, and quite often one encounters particular results for the equilibrium distribution of a continuum of charges scattered in the literature.
The aim of this work is to bring together the existing results so that they can be easily consulted.In order to pursue this task we will proceed in a constructive fashion: starting from known results in the theory of elasticity, we will rigorously prove an exact one-to-one correspondence between the physics of charges in equilibrium and the formalism of contact mechanics.During this endeavor, we shall also provide the required formalism for dealing with problems that may require numerical analysis.In this way, known results in contact mechanics will automatically become what one requires for the study of equilibrium charges under different types of confinement.Again, this connection is not new, but since it is so deeply rooted in the jargon of potential theory, it is basically unknown to most authors.
Additionally, the particular case of particles bound to move on axially symmetric surfaces will lead to new and interesting results regarding the ensuing equilibrium charge density distributions.In addition, the modification of the Coulomb interaction will also be shown to provide in some cases analytical solutions that generalized previous particular instances.
The present contribution is divided as follows.In the first section we shall discuss the variational formalism for finding the equilibrium density of interacting charges under a general confinement.Next, the so-called Boussinesq problem [32] in contact mechanics will establish a parallelism drawn from the aforementioned field of physics and electrostatics, that will allow us to establish a bijection between equilibrium elastic mechanics and the equilibrium distribution of interacting charges.The most common cases of confining potentials will then be presented, as well as the reference to others.The following sections will be devoted to the exploitation of the variational formalism in order to find equilibrium charge density distributions in several axisymmetric problems with a hard-wall confinement.We will provide new numerical approaches for the finite case, as well as the solution of integral equations for the continuum extension.The case of a power-law extension of the Coulomb interaction will allow us to find analytic results in the generalized d-dimensional ball.Finally, some conclusions will be drawn.

Solution of Equilibrium Problems in Electrostatics
A successful approach to the problem of finding the equilibrium configuration for a finite number of particles N ≫ 1 requires considering the system in the continuum limit as a first approximation.Further contributions account for finite-size effects, which are sometimes called "correlation energies."They are mostly estimated via local density approximation techniques, [9] that assume that for clusters of large enough charges the density can be understood to be constant.
The usual procedure is to solve either the Laplace or the Poisson equation in the appropriate coordinates for the potential V(x, y, z), and then retrieve the (usually surface) density as that is, the negative normal derivative of the potential.Now, the complexity of the problem then lies in solving the corresponding Laplace or Poisson equation analytically or, in the worst case, to numerically solve it by means of finite difference methods.Another approach to solving the same problem relies on applying the variational principle to the total electrostatic or self-energy of the system in order to obtain the unknown charge density.Returning to the previous example of the disk, we want to obtain the unknown charge density  H (r).
For charges confined by a hard-wall potential, the energy expression in Equation ( 1) can be approximated in the continuum by the integrals over the disc The continuum approximation treats the density  H (r) as a smooth function rather than the sum of delta functions where r i is the position of the ith charge.One then minimizes the energy of the cluster, with respect to the smooth function  H (r), subject to the constraint that the number of particles is constant.Introducing the Lagrange multiplier  the constrained equation is A variation in the energy is given by where  H (r) represents a small change in the charge density.
Keeping only terms up to first order, Equation (8) gives To make the functional derivative stationary, we require Now, we have arrived at an integral equation for the density, depending on a Lagrange multiplier.To solve this integral equation it is convenient to write the integral in Equation (10) in terms of radial and angular variables.
where K(k) is the elliptical integral of the first kind.This is a Fredholm integral equation of the first kind, which shall possess a radially symmetric solution.Admittedly, solving an integral equation is not a straightforward endeavor, nor to approach the solution of particular potentials V(x, y, z).The closed solution of the previous integral equation (IE) can be obtained by the method of transforming operators and expressed by two quadratures, and it is illustrative to present the details as compared to the complexity of usual electrostatic potential-solving methods.Based on the relation between the Abel rotation operators and the operator generated by the singular Cauchy kernel, solving this IE can be reduced to solving the simplest singular integral equation, again admitting a closed solution.Using the well-known spectral relations containing Legendre polynomials, its analytic solution can also be expressed in the form of a series in these polynomials.From the point of view of the numerical analysis of the problem, the method of Legendre orthogonal polynomials seems to be quite simple and effective.In this method, the solution of the IE for an arbitrary continuous function (r ′ ) is shown.Also, it clearly shows how to proceed in an eventual numerical computation.In such a case, one can gain sufficiently high accuracy by reducing the concomitant IE with the subsequent application of a wellknown numerical-analytical method based on Gauss quadrature formulas for calculating integrals at Chebyshev nodes.
After introducing the new set of variables r = R,r ′ = R (0 < ,  < 1) the previous IE reads as subject to the condition Next, we write down the well-known spectral relations [33,34] ∫ where P 2n (x) are the Legendre polynomials of index 2n, and Γ 2 (n+1) , where Γ(x) is the known gamma function. [35]We also write the orthogonality condition for the Legendre polynomials with the new variable x = √ 1 −  2 as Now, proceeding from (15), we present the solution of IE (13) in the form of the infinite series with unknown coefficients x n .After some algebra, we obtain After using the previous orthogonality relations and extending f () evenly to the interval (−1, 0), we find These coefficients are easily calculated using the Gauss quadrature formula at the Chebyshev nodes.After some substitutions, x 0 = f 0 ∕ 0 and therefore which ultimately establishes the relationship between R and N.
In our particular case when we have and consequently Incidentally, this result is the density profile obtained if a hemispherical shell of charge is projected onto a plane, a result known to Green, Kelvin, and Weber. [36]Usual known equations for equilibrium densities are obtained in the literature as limiting cases of the ellipsoidal shell as an elliptical hard-wall confinement.Now, suppose that a general potential V(r) is applied to the system, again in the setting of the 2D planar disk.There is no hard-wall confinement now.V(r) is responsible for holding the particles together against the Coulomb repulsion.The total energy functional now reads Proceeding as before, the ensuing integral equation reads as which is equivalent to (10) when an external potential is applied.Now, rearranging (25), we have which is nothing but a well-known result in potential theory.Therefore, if we are able to solve (26) for an arbitrary potential V(r), we have achieved our goal.

The Boussinesq Problem in Contact Mechanics
The Boussinesq problem [32] happens to be known as the situation when one is interested in knowing the pressure distribution As defined in Podio-Guidugli and Favata, [37] "the Boussinesq problem (Joseph Valentin B., 1842-1929) consists in finding the elastic state in a linearly elastic isotropic half-space, subject to a concentrated load applied in a point of its boundary plane and perpendicular to it."An indenter on an half-infinite elastic media exerts stresses normally (pressure − zz and perpendicularly  zr ) to the point of action.The mathematical dependence on r of the action at a distance r perpendicular to the indentation is what makes the problem relevant mathematically.See text for details.
when the surface of an object acts normally (indenter) on the semi-infinite surface of an elastic media.Originally in its simple form, it corresponded to the computation of the reaction in the elastic media of the action of a point-like pressure P, as shown in Figure 1.What is relevant, regardless of the details of the problem, is that the solutions of the normal displacement constitute harmonic functions, and as such, satisfy the Laplace equation.Thus, when the normal surface (z = 0) displacement is obtained as with E being the Young modulus and  the Poisson number, there appears to an obvious similarity between (27) and the electrostatic potential created by a point charge.This fact precisely constituted the leitmotif for extending the potential theory developed in gravitation to the both contact mechanics and electrostatics.We have not derived or used the list of equations leading to the final expression (27).The rationale is that we thus obtain the main ingredient for the forthcoming bijection, namely, the 1∕r dependency of an infinitesimal pressure distribution on an elastic media.Thus, it is plain that the sum of the contributions of all normal displacements (i.e., minus the component of the stress  zz corresponding to an infinitesimal of the pressure distribution under the indenter) shall be mathematically equivalent to the contribution of the total energy arising from the Coulomb interaction (potential) between infinitesimal charges.
The Boussinesq problem is the conventional nomenclature to introduce the effect (normal deformation) that produces a concentrated normal force acting on an elastic half-space.But owing to historical rigor, we should have the following timeline: [38] i) Coulomb discovers the 1∕r dependency of the electrostatic potential by the time of the French Revolution; ii) Poisson generalizes the Laplace equation to arbitrary mass distributions, extending Laplace's methods on potential theory; iii) almost simul-taneously, Poisson worked on equivalent problems of charge distributions, giving birth to modern electrostatics.iv) Although Poisson failed to solve the potential for the charged ellipsoid, Prana did in 1840 and was later employed by Hertz in his seminal contribution to the contact problem between two elastic bodies, therefore founding the field of contact mechanics.Incidentally, the results of Hertz predate those of Boussinesq by a few years.In any case, from a historical point of view, the merit of Hertz's noticing the similarity between electrostatic and elastic problems is precisely what we revisit in the present contribution.
The classical potential theory apparently goes back at least to Gauss, Green, Weber, and Kelvin; other notable contributors (in roughly chronological order) include Poincaré, de la Vallée-Poussin, the Riesz brothers, Leja, Frostman, Szegö, Keldysh, Lavrentiev, Brelot, Ohtsuka, Choquet, Fuglede, Rakhmanov, Totik, Saff, and Dragnev.The field has a number of important connections to other areas of mathematics, for example, theory of PDEs, polynomial approximation, martingales, theory of holomorphic functions, measure theory, and so on.
Posterior work [38] mainly consisted in developing techniques for solving pressure/charge distributions analytically and numerically, up to a point where specific results in contact mechanics had no counterpart in equilibrium electrostatics (theory of cracks, adhesion, etc).
The continuous transfer of knowledge from one field of physics to another can be seen purely on mathematical grounds from potential theory, but this approach lacks the perspective than most physicists have due to the degree of specialization of the different branches of physics.Therefore, highlighting the known results in contact mechanics is of great utility in equilibrium electrostatics, which shall be explored in the following.

Bijection between the Axisymmetric Problems of Electrostatics and Contact Problems of the Theory of Elasticity
One of the most important cases will correspond to planar problems where the topology of the system displays a symmetry around an axis, although the bijection extends also to non-axially symmetric cases due to the same mathematical structure underneath contact mechanics and equilibrium electrostatics, namely, classical potential theory.Incidentally, another area of application of potential theory is found in heat transportation, when the flow of heat is in equilibrium.
The boundary value problems of electrostatics and the classical axisymmetric (and non-symmetric) contact problems of the theory of elasticity, based of the hypotheses of Hertz, [39] are closely related by the general ideas of potential theory and the mathematical equations that describe them.Therefore, under certain conditions, it is possible to establish a one-to-one correspondence between them.For this purpose, we briefly outline the formulation of the classical contact problems of the theory of elasticity and derive their governing IE.Following Shtaerman, [40] we will consider the problem of compression of elastic bodies, initially touching at one point.The origin of coordinates O of the right rectangular coordinate system xyz is placed at the initial point of 0 of tangency of the compressible bodies.Let z = f 1 (x, y) and z 2 = −f 2 (x, y) be their surfaces, which we assume to be smooth enough to be continuously differentiable sufficient number of times.In addition, we will assume that the elastic bodies are compressed against each other through vertical concentrated forces of magnitude P, such that the bodies can only translate in the vertical direction.Then, in the process of compression, a contact area is formed around the initial point of contact , the geometric shape of which, generally speaking, is unknown in advance.
For a mathematical description of the stated problem, we use the Hertz hypotheses, [39] namely, i) the dimensions of the contact area (pad)  are much smaller than the typical dimensions of the compressible bodies; ii) due to the smoothness of the surfaces of the bodies in contact, only normal stresses act on the contact area; and iii) the elastic displacements of the compressible bodies are negligible as compared to the dimensions of the contact area.By virtue of hypothesis iii) for vertical elastic displacements  1 and  2 on each corresponding body, the following condition must hold which is called the contact condition.Here  is the mutual vertical convergence of the bodies during their translational movement, which characterizes the rigidity of the compressible bodies and is subject to determination.Further, we denote by Σ the projection of the contact area  onto the horizontal plane xOy.According to hypothesis ii) at the point (x, y) of this area Σ only normal pressure p(x, y) applies.Now, by hypothesis i) we replace the compressible elastic bodies with upper ( 1 =  1 (x, y)) and lower ( 2 =  2 (x, y)) elastic half-spaces caused by the contact pressure p(x, y) in the region Σ located at z = 0. Using the solution of the well-known Boussinesq problem [41] and the linear principle of superposition in the theory of elasticity, we will have (E 1 is the elastic modulus and  1 is the Poisson ratio of the upper half-space of the first elastic body; likewise (E 2 ,  2 ) for the second body).The previous contact condition will now read as from which the contact pressure p(x, y) is determined.The solution to IE (30) must satisfy the condition which implies equilibrium of each of the two compressible bodies.Usually, the first body is considered absolutely rigid, which implies IE (32) and condition (31) constitute the link between contact mechanics and the classical potential theory.Namely, we introduce into consideration the potential of a simple layer with the density of sources p(x, y) on a double coated flat surface Σ.The potential u(x, y, z) outside the flat region Σ in the entire space satisfies the Laplace equation, and at infinity, according to (31), has the behavior Consequently, the solution of IE (32) under condition ( 31) is equivalent to the solution of the following regular exterior Dirichlet boundary value problem After solving problem (35), the density of sources p(x, y) will be determined by the well-known formula for the normal derivative of the potential of a simple layer which nicely coincides with the equivalent (3) case in electrostatics.Figure 2 describes the correspondence between contact mechanics and equilibrium electrostatics.Thus, we have reached a point where we can finally establish a bijection between the results in equilibrium electrostatics and contact mechanics.If the axisymmetric case is finally given by the integral expression then the governing IEs of the electrostatic problem and the corresponding condition shall read as Since (r) is a Lagrange multiplier, then one must have (r) =  = constant and, therefore, the right-hand side of IE ( 38) is represented in the form These problems are physically united by the general ideas of the theory of classical potential, and mathematically described by the same equations.Namely, • the density of charges in the general case (r) corresponds to the contact pressure p(r) under the indenter, and vice versa; • the total value of the charges N corresponds to the vertical force P pressing the indenter to the elastic base and vice versa; • the confining potential V(r) corresponds to the function f (r), which, up to a constant, geometrically describes the base of the indenter.
This correspondence (as depicted in Figure 2), which is known to occur for all harmonic functions under the similar Dirichlet problem (35), has benefited both fields of study, quite often complementing each other as far as different shape indenters (or confining potentials) and underlying assumptions are concerned.We thus bring this correspondence to the general audience so that it is unnecessary to solve integral equations from scratch.
Once the bijection is clear, we shall report the results corresponding to the pressure distribution p(r) in contact mechanics as the electrostatic equilibrium surface density (r), corresponding to different confinements V(r) and V(r, ).In the following, we review the solution for several confining potentials.There is an extensive literature that discusses these results.We refer to refs.[40, 42-49] and references therein for instances not covered here.The main source of references for the following confinements is Popov et al. [42]

Linear Confinement
We have already obtained the result for the equilibrium surface density inside an infinite circular confinement.There is a nice physically motivated case of a radially applied electric field in the circular hard-wall confinement.In this case, it is known that the radial symmetry is preserved, V(r) = |qF|r (from the equivalence with contact problems).Thus, the potential V(r) is equivalent to that of a linear indenter up to a certain radius R and, therefore, the ensuing equilibrium distribution is of the form which clearly diverges at the origin.

Parabolic Confinement
The application of a magnetic field perpendicular to a circular disk and/or having it to rotate, induces the combined potential V(r) = (B + m e e )  2 r 2 .Notice that the null magnetic field case contributes very little, unless  is large.The original solution to the parabolic confinement goes back to the classical work of Hertz, [39] who used potential theory, and reads as again up to a certain radius R.

Power-Law Confinement
The general case of the type V(r) = c r n up to a certain radius R, with n being a positive real number, is given by where B(⋅; ⋅, ⋅) denotes the incomplete beta function B(z; a, b) ≡ z a a 2 F 1 (a, 1 − b; a + 1; z), with the latter being the hypergeometric function.

Infinite Elliptical Confinement
This instance was already known to Green, Kelvin, and Weber, as it can be obtained as the projection of a thin ellipsoid on the z = 0 plane.It is given by Here, a and b are the typical semi-axes.

Confinement that Induces a Constant Charge Density
It is possible to design a confining potential in such a way that the induced charge density is uniform.This problem was initially solved in terms of contact mechanics, but using potential theory, by Boussinesq and extended by Lamb [50] in 1902 in the form of hypergeometric functions.The following derivation is due to Föppl. [51]Again for a certain radius R, we have Here E(⋅) is the usual complete elliptic integral of the second kind.

Infinite Annular Confinement
The closest instance in equilibrium electrostatics corresponds to a hole in an infinite planar conductive media, [52] which occurs in an asymptotic limit of the outer radius b → ∞.This is well-known solution that can obtained by a recourse to standard electrostatics.However, the instance where both the inner a and outer b radii are finite does not admit a closed form solution.Arguably, the general strategy in treating problems with annular regions is to extend the "annular" condition either to r = 0 or to r = ∞ and relax the consequent overprescription by superposition of an additional two-part solution.
The ensuing Fredholm integral equation is solved iteratively and is given as power series in terms of the ratio b∕a, as given in ref. [53].In the limit b∕a → 0 one nicely recovers the charged disk in equilibrium.A suitable intermediate approximation is usually given by We can extend the b → ∞ case when it is wrapped around itself, that is, a small spherical charged cap on the S 2 sphere.There is one way to address the problem for small aperture angles  c in spherical coordinates.The idea is to project or extend the planar solution to the sphere for a finite curvature.Let us call R ′ the maximum distance between the pole of the cap and any point at  =  c .Basic trigonometry implies that R ′ ≡ 2R sin(  c 2 ).As expected, when the radius of curvature tends to infinity,  c tends to zero so that R ′ remains constant.In other words, lim R→∞ R ′ = lim →0 R ′ = R || , where R || stands for the usual radius of the circle in the planar case.Since planar distribution on the disk goes as 2 ) and the variable r by the distance d ≡ 2R sin(  2 ).Thus, we have (R Therefore, after imposing charge normalization, we obtain that the equilibrium surface charge density for small angles {,  c } on a spherical cap goes as One can easily check the validity of ( 46) by computing the total ensuing self-energy and observe that it remains lower as compared to the uniform density case for a broad range of angles (they coincide around the hemisphere case, that is, when  c = ∕2).

Null Inner Potential Followed by a Linear or Quadratic Contribution
In the language of contact mechanics, we are referring to an indenter with a flattened tip followed by a i) linear profile (truncated cone) or ii) a parabolic profile.A confining potential V(r) is null up to a radius b, and then changes linearly or quadratically (truncated paraboloid) up to the radius a.This combination is certainly interesting for it has been considered in contact problems only.We have [42] • "Truncated cone" • "Truncated paraboloid" As explained in ref. [42], the former expressions can in principle be expressed in closed form by using elliptic integrals, as in the original computations.However, the resulting expressions are very extensive and cumbersome to evaluate, which is why a numerical evaluation will definitely be preferable.

Undulated Potential
It is possible to induce an equilibrium oscillatory charge density profile, an electrostatic analog of optical trapping, along the radius of the confinement.The earliest reference to such problem is due to Gallop [54] (1886) and MacDonald [55] (1895), although their original results where simplified by Copson [36] in 1945, all of them working in electrostatics and potential theory.The confining potential is of the type V(r) = J n (cr), that is, the Bessel functions of the first kind with n being an arbitrary positive integer and c a positive real constant, which modifies the scale of undulation.
In this case, the general solution is provided by Thus, one gets, besides the original confining potential, an additional oscillatory term.A more recent approach was made by Guduru [56] in a context of contact mechanics.Given a confining potential, with h being the amplitude of the oscillation and  its wavelength.The ensuing equilibrium density distribution is given by Here R is the maximum radius until which the potential acts.
H n (⋅) denotes the nth order Struve function, which is related to the generalized hypergeometric function 1 F 2 by ).

Polar Potential
Let us apply an external electric field F along the x-direction of a circular disk.The contribution as in (25) will be V(r ′ , ) = ±|e|Fr ′ cos , depending on the sign of the charges, but in the end the equilibrium density will possess the same symmetry around  = 0. Thus, we must seek the solution to the following integral equation over the disk domain Ω.It is apparent that (r) = (r, ).The solution is known [36] to be Also of interest is the case of an ion trap, where electric quadrupoles can generate polar potentials.Once the disk is placed in the center of coordinates, we can apply two types of in-plane confinements.The first one is V(r, ) = A(x 2 − y 2 ) = Ar 2 cos(2), with solution [36] (r, ) = 8A r 2 cos(2) The second quadrupolar case corresponds to V(r) = A(x 2 + y 2 − 2z 2 ) = A(r 2 − 2z 2 ).Since z = C is a constant, only the radial part intervenes, and we recover the ∝ r 2 confinement.

Infinite Potential with Non-Circular Boundaries
All polygons can in principle be solved analytically, usually by means of a power series plus a contribution that takes into account the divergence of the equilibrium distribution on the boundaries.However, all of these instances share the same problem, namely, the concentration of charges at the corners.This problem was first addressed by Noble [57] in 1959.He used spherical coordinates and was led to an eigenvalue problem which involved mixed boundary conditions, since the plane angular sector corresponds to only part of a coordinate surface.Noble basically obtained approximate solutions to the resulting dual equations by means of an integral equation approach with the introduction of an eigenvalue problem that is solved by depending on the angle sector.This method was later extended by Morrison and Lewis [58] (1976).Different cases are the ones where the potential V(x, y) is a linear or polynomial function of x and y, which are solved for any polygonal shape.There, the results are analytic.An illustrative case is provided in ref. [59] and the references therein.One has to bear in mind that the former cases have their origin exclusively in contact mechanics.Important instances with hard-wall confinements constitute the triangle, the rectangle, and the hexagon.The case of the rectangle, first studied by Love [60] is semi-analytic, accounting for the divergence at the boundaries but only admitting an infinite power series.The triangle is given in similar terms by Svec and Gladwell. [61]In general, all distributions diverge at the boundary Ω, even if the pressure is uniform, as well as their derivatives, as some fractional power law (Figure 3).
We shall consider a direct minimization of the total energy for several boundaries of charges confined in hard-wall regions.The accumulation of the contours described by their equilibrium configurations provides a sense of how the true equilibrium distributions behave.

Other Instances Related to Equilibrium Electrostatics Only
The present bijection between contact mechanics and equilibrium electrostatics allows us to utilize results well known in the former for the latter cases.Thus, the common background of potential theory is avoided.Now, the same tools introduced before (mainly the variational principle) can be easily extended to problems in equilibrium electrostatics that are not found in the literature.As we shall see, these instances serve the purpose not only of a natural arena where to apply and extend previous methods, but to round up the current work with novel results.
Let us consider other hard-wall problems, only differing in the topology of confinement.There is also a hard-wall confinement, but the regions where the charges are deposited will vary.That is, we shall have N interacting charges with total energy where V(r), the potential energy, is zero in the domain Ω and infinite on the boundary Ω.The combined nature of both Ω and Ω will define the cases under scrutiny.
In order to pursue the equilibrium distribution of the continuum (r), we shall choose an axially symmetric confining surface given by z = f (x, y) = f (r).The surface element is then given by The equivalent of the general integral Equation (10) shall read as Since we shall be able to set r 1 = 0, the upper limit in (55) will be given by a typical radius R. Equivalently, the normalization constraint reads as Once we have the general setting, we can address several cases, such the following ones: • the cone.For z = a r, we have h(r) = √ 1 + a 2 • the circular paraboloid.For z = a r 2 , we have h(r) = √ 1 + 4a 2 r 2 • the general surface of revolution z = a r n • the circular ellipsoid (cap).For z = b The aforementioned cases lead to integral equations, whose solutions for the density distribution (r) are the ones we are looking for.A particular case of interest is given by the cylindrical shell.The equivalent of the general integral Equation (10) for the present problem reads as where K(k) is the elliptical integral of the first kind.We have made use of the fact the the surface density depends only on the variable z due to symmetry arguments.It is apparent that the Coulomb interaction appears convoluted by the action of several rings of charge.The solution of Equation ( 57) is subject to the constraint (6) in the form Due to the symmetry of the problem, we must expect (z) to be symmetric around z = L∕2, that is, the midpoint between both ends of the hollow cylinder.A zeroth approach to obtain an analytic solution of (57) is to consider the limit of very large L, that is, |z − z ′ | ≫ R.Then, the integral equation becomes nothing but a convolution of the density (z).Then, one Laplace transforms both sides of the approximated IE and finds (z) after a non-trivial inverse.However, the result is analytic and given by a uniform density (z) =  u .

Numerical Approach to Finite Number of Charges: Three Methods
There are three approximations to numerically retrieve an equilibrium density distribution of N charges.The first method consists in a direct minimization via heuristic algorithms, and is relevant if one is interested in studying finite-size effects, such as the formation of clusters.This method is straightforward but requires considerable computational power.The second approach relies of the symmetry of the problem.It is still finite, but in those cases where we have rotational symmetry, one can introduce rings of charges, alongside with some shell-occupation numbers.In each case, a numerical optimization of the total interaction energy was obtained by an interior point optimization algorithm [62] with inequality constraints.When the number of charges N ≫ 1, one can observe a flat and compact distribution in the innermost regions, which implies an almost uniform charge density away from the boundaries.Close to the boundary, on the other hand, accumulation of charges indicates a spike in density.See text for details.The interaction between rings is averaged over all possible angle offsets.Finally, one has the option of solving an integral equation numerically.Let us explore in the following the aforementioned techniques.
One numerical approach to the equilibrium configuration of N charges distributed will be that of the surface cylinder with open ends.A simulated annealing Monte Carlo method has been performed for N = 600 in Figure 4 for R = L = 1.Notice the formation of rings.In such a case, one can take advantage of the symmetry at both sides of z = L∕2 combining the numerical results with a restriction on the positions of the shells z i , thus incorporating z i as optimization variables and greatly reducing the overall number of computations.
For N = 600 charges, we have obtained a total of 10 rings for d∕R = 1 (for instance, we have 7 rings for N = 300), as shown in Figure 4.The occupation numbers are {110, 58, 48, 41, 41, 42, 44, 50, 57, 109}.Then numerical equilibrium configuration has slight localized ripples, which is typical for high-N optimization.The ensuing total energy is 182 513.734.Once the ripples are smoothed out, the energy becomes 182 517.701.
Although the equilibrium density we are seeking is the one corresponding to the continuum, in the discrete case the charges will distribute themselves on rings, where each ring will contain its particles evenly distributed on the periphery.As mentioned, all these rings will be distributed symmetrically around the middle of the length of the cylinder.Thus, one can assume from the very beginning that these structures will be the ones occurring in practice.The goal is to avoid extremely long optimizations of the total energy for each charge.In this fashion, one can approximate the equilibrium density distribution (z) obtained by numerical means.
The second approach will arise from introducing the selfenergy of each ring plus the contributions from all interactions between rings, averaged.The present approach is a modified one from the original idea of Cerkaski et al. [63] The number of charges N shall distribute themselves into s shells or rings, and the corresponding occupation number shall be denoted by n.This problem is reminiscent of the one in mathematics corresponding to a integer partition of a number N into a sum of s terms, that is, N = n 1 + n 2 + ⋯ + n s .The problem here will be to seek several partitions with an unknown number of shells or terms s, and to choose the one that minimizes the total energy.Also, each occupation number shall correspond to a distance z i along the axis of the cylinder.In the case of Figure 4, when we approach the problem by averaging the energy with regard to relative angular positions, that is, with unknown positions of rings z i and occupations n i , we obtain an averaged energy of 182 522.023.The concomitant n i and z i ∕R are {110, 0.; 57, 0.0930476501; 49, 0.203985325; 41, 0.31915729; 41, 0.434227197; 42, 0.555223433; 44, 0.676873087; 49, 0.797795154; 57, 0.908855614; 110; 1.}.All in all, the occupation number approach retrieves good approximate results for high N, with fewer unknown variables.
Let us consider the general case of a surface given by the revolution of the curve z = f (r) around the z-axis.The total (averaged) energy can be shown to read as with Δz ij ≡ |f (r i ) − f (r j )|.Distances could be rescaled in terms of some characteristic length L. The nature of each object shall be encoded in the relation z = f (r).
In the specific case of the cylinder, with Δz ij = 0, the total (averaged) energy of the N charges distributed in the previous setting shall read, in units of k e e 2 ∕R, as Positions along the axis of the cylinder z ′ i,j are given in units of the radius R. K(k) is the usual elliptic integral of the first kind.The case of the cylinder will be the only one with two open ends.Now, a complementary condition that may help retrieve the positions of the shells is provided by where the positions z are given in units of the radius R of the cylinder, and E(k) is the elliptic integral of the second kind.The set of nonlinear equations thus defined is to be solved for z k and the occupation numbers for each shell n k .The last numerical method involves the solution of an integral equation.Let us thus consider the problem of finding the surface charge density in equilibrium for a cylindrical shell.The extended problem of finding the charge densities when a conducting cylinder is closed with its bases has been the subject of intensive work.As Smythe remarked in his pioneer work, [64] "the literature is blank on this subject."He addressed the problem by introducing two densities, namely, one on the bases and one for the cylindrical shell.These densities complied with the symmetries of the problem and imposed an ad hoc restriction of their functional form (special divergence of the kind  − 1 3 where  is the distance from the edge, identical to the one expected from a rectangular wedge).Despite the efforts, it still remains an open question to find a closed form for those equilibrium densities.
Since we are going to address a cylindrical shell only, we expect a divergence of the  − 1 2 kind.Let us consider the original integral Equation (57) subjected to the restriction (58) in the following new variables and the condition To the best of our recollection, (62) possesses a priori no analytic solution.Thus, we proceed with the numerical procedure.Let us introduce the difference kernel As  → , k → 1 and the complete elliptic integral of the first kind K(k) diverges logarithmically. [35]Proceeding from this, we reduce IE (62) under the condition (63) to a singular IE, separating the principal singular part of the kernel.Specific singular IE shall follow other kinds of prescriptions.In our case, we shall proceed to differentiate the function K(k) and calculate where E(k) is the complete elliptic integral of the second kind of module k.By rearranging the previous expression, we have What one achieves is to split the L 1 ( − ) kernel as the sum of its principal singular part in the form of the so-called Cauchy kernel with coefficient −1∕2, and its regular part in the form of the remaining expression.By taking this circumstance into account, the IE (62) will now read as the following singular IE: under condition (63).The previous procedure ensures that the kernel M( − ) is bounded at 0 < ,  < .For our purposes, we shall extend the integral intervals to the standard (−1, 1) by setting Finally, our singular IE problem will take the form Now, the standard procedure in these cases [65,66] is to represent the solution of ( 69) by the formula where Φ(u) is the Hölder function on the segment [−1, 1].We now choose an arbitrary M and determine the Chebyshev nodes Here u m are the roots of the Chebyshev polynomials of the first kind T M (u), while t r are the roots of the Chebyshev polynomials of the second kind U M−1 (t).
The ensuing IE problem is now reduced to solving a system of linear algebraic equations After solving (72), the solution at the nodal points is computed according to Higher number of equations M provide a better precision for the function (u) or the original sought function (z).
We have applied the previous numerical recipe to the cylindrical case for several values of  ≡ L∕R and a system of M = 1000 linear equations, which is solved for arbitrary .The numerical solutions of the corresponding integral equation match perfectly the functional form (u) = 2  (1 − u 2 ) 1∕2 .We not only find the solution to the equilibrium charge density distributions, but in the case of the cylinder with open ends it is found to be analytic.Remarkably, lim →∞  ⋅ (u) = 2∕, which nicely recovers the known result when L is large as compared to R. Conversely, (u) is not well defined as  → 0 + , for the total energy diverges logarithmically (we recover a charged ring).When the density is considered to be uniform, which constitutes an upper bound to the actual equilibrium energy, it is found [67] in an analytic fash-ion: with  ≡ L∕R.
For the sake of completeness, we have also solved the integral equation corresponding to the cone with no base, for several values of H∕R, as well as for the hemisphere (circular ellipsoid with a = b = R).The results are depicted in the inset of Figure 5.The minimum of the density shifts from left (H∕R = 1∕2) to right (H∕R = 5) in steps of 1/2.Notice the competition between the accumulation of charges at the vertex of the cone with those at the periphery of the base.In the limiting case of H → 0 + , one recovers the circular disk equilibrium density ∝ [1 − (r∕R) 2 ] − 1 2 , as can be appreciated from Figure 5. Regarding the hemisphere, the equilibrium density as measured perpendicularly to the North pole, is mostly uniform until is diverges (numerically checked in logarithmic fashion) as r → R.
It is thus plain that the equilibrium surface charge density of any axially symmetric object with a hard-wall confinement, that is, with particles bound to move of the surface, can be obtained via the numerical or analytic solution of a Fredholm integral equation of the first kind.In most cases, the ensuing IE will be singular.Thus, according to its nature, it shall be solved accordingly.

A Novel Approach: The Disk with Two Charged Species
A simple calculation with uniform charge distribution shows that, for equal densities, the particles with more charge tend to concentrate on the periphery, whereas the "lighter" ones are pushed inward.Remarkably, for a finite instance of N∕2 particles each, the finite-size effects are very much present.Only in the limit of N ≫ 1 one is expected to reproduce two clear regions.
There is very little literature on the subject. [68,69]Barber [68] describes the problem by an annular discontinuous punch, which is mathematically equivalent to having an electrostatic potential field due to a circular disc and a surrounding annular disc, maintained at different potentials.But this is not exactly our problem.Barber [68] proceeds to solve the problem containing two simultaneous Fredholm equations for two unknown functions, which has to be carried out numerically.In any case, the region fixes the separation between the inner circle and the annular region, while it should depend in the discrete case (and the continuum case as an extension) of the ratio of charges.
Thus, in this case we are left alone with the numerical results for the discrete system (Figure 6).

Important Instance: The Circular Capacitor
Suppose that we have two charged disks with opposite charges.The first disk lies on the xy plane with charge −Q, and the second one a distance d above being parallel and coaxial, with charge +Q.Due to symmetry reasons, the functional form for the distribution of charges in both disks is the same, only differing with Figure 6.Distribution of N charges in equilibrium for two species of charges q, blue, (N∕2) charges and 2q, red (N∕2 charges), for N = 20, 50, and 1000.For N small, one observed an equilibrium that mixes the positions of both species.As N grows, there appears a clear distinctive pattern.In the continuum, it is energetically favorable to have the two species spatially split.The ensuing densities diverge additionally at the critical radius R. We are unaware of an equilibrium problem of this kind in electrostatics that mixes two differently charged species.Thus, it is expected that R depends solely on the quotient of charge species.In the rightmost case, it is roughly R ≈ 0.8.See text for details. the sign.Since we are looking for the equilibrium charge distribution, we can safely assume that the uniform charge density distribution will constitute a reasonable upper bound to the real minimum energy value.
A direct (discrete) treatment of the problem is given on the left in Figure 7, where the minimal configuration for the total energy of 1000 charges and a distance d∕R = 0.5 between plates is shown.It is contrasted with the minimal configuration produced by the single-species energy for s = 1 and N = 1000 on the disk.The middle graph shows smoothed density histograms for the two configurations, in blue for the optimal capacitor distribution, and in red for the single-species distribution.It can be seen that the single-species energy leads to higher density toward the boundary of the disk.This phenomenon is in agreement with the capacitor case corresponding to the less strongly repulsive kernel, see Equation (76).
In the continuum case, the total energy (up to an additive term if a constant potential V is applied between plates) shall read as Taking into account the equivalence between the functional forms for the densities, the energy becomes Then densities (r) from both plates are subject to the constraint that the total number of particles is constant and the system is electrically neutral (|e|N∕2 − |e|N∕2 = Q − Q = 0).Proceeding as before, the functional derivative being stationary requires Writing the previous integral equation in terms of radial and angular variables we have The integral Equation ( 79) is to be interpreted as the total potential on the upper disc due to the contribution from the positively charged (upper) disk plus the one from the negatively charged disk (lower one).Notice how the limit d → ∞ recovers the problem for one single disk, whereas for d → 0 the problem disappears for the system becomes electrically neutral.
We have solved the integral Equation ( 79) for the equilibrium charge density (r) and several distances d between plates, and the solutions are depicted in Figure 8.
The equilibrium distribution tends to the well-known case of the disk as d → ∞, which numerically occurs around d∕R ≈ 10.
[72] The second term in ( 79) is related to the potential between plates, responsible for the same amount of charge Q deposited in each disk, but with opposite signs.In principle, integral Equation ( 79) could be solved analytically if we allow a series expansion in terms of d.Thus, it is straightforward to compare with the most recent progress in this regard.
As far as the numerical treatment is concerned, a simple numerical exploration shows that the minimum energy for a system of positive charges on the upper plate, plus the system of negative charges on the lower, plus the interaction between them is achieved when the positions on both disks are mirrored.This fact implies that the rings formed by the charges, as well as their occupation numbers, will be exactly the same on both plates.The total energy in the ring approach with average interaction up-up, down-down, and up-down will thus become The first k = 0 term is responsible for the Coulomb interaction between particles only differing by a distance d (one on top of the other).
The usual limits are recovered: for d → ∞, one is left with twice the averaged energy for a single disk.On the other hand, for d → 0 one has to be careful: whereas the inter-ring energy disappears, a term − N∕2 d is left (the aforementioned k = 0 term).This term is present because opposite charged particles are vertically arranged in pairs, but in the limit of null distance this is most certainly not the minimal configuration.There, opposite charges would occupy diametrically opposed positions and their contribution would become identically zero.The expected equilibrium (r) shall not greatly differ from the one for a single disk.

Stored Coulomb Self-Energy of Some Uniformly Charged Bodies
Retrieving the equilibrium charge distribution is not always possible in an analytic fashion.Even when it is so, the computation of the total self-energy of the body may entail a numeric calculation after all.Thus, in some cases when certain dimension of the body is large, one can assume the charge density distribution to be uniform.Doing so greatly simplifies the problem and, in some instances, the total Coulomb self-energy can be computed analytically.What is of importance in the sort of cases is that the stored energy for the uniform case always acts as an upper bound to the real equilibrium one obtained by the methods of previous sections.In the following we show the energies of previous instances and comment on the results as compared to their equilibrium counterparts.
• Uniformly charged rectangular plate.The exact calculation of U(L x , L y ) is very challenging with the final result obtained by Ciftja [73] that reads: • Uniformly charged equilateral triangle.The final expression for electrostatic self-energy of a uniformly charged equilateral triangle is [74] E u = 2 ln( 3) where L is the length of the equilateral triangle.• Uniformly charged annulus.The final result for the selfenergy is given in the following form: [75] U • Uniformly charged cylindrical shell.Notice that the total energy of a uniformly charged cylindrical shell is known in the literature, and that it decreases rapidly with  (being divergent at  = 0).This example shows how different the two instances are.The final result for the Coulomb self-energy is [67] U with  ≡ L∕R.• Uniformly charged hemispherical surface.As shown in earlier work, [76] the final result for the self-energy is remarkably simple:

The Circular Case with 1∕r s Repulsion between Charges
The case of a circular hard-wall confinement with 1∕r s repulsion is the closest analog of the solid sphere treated by Mughal in ref. [77].The ensuing integral equation is with s being a positive value.According to Hardin and Saff, [78] there are two regimes to be considered, namely, i) 0 < s < p(= 2) (not uniform), and ii) s ≥ 2 (uniform).The so-called "poppy seed bagel theorem" (Hardin and Saff [78] ), a name that alludes to discrete equilibrium configurations on the torus, according to which the configuration of points that minimizes the Riesz energy ∑ 1≤i<j≤N 1∕|r i − r j | s on a rectifiable manifold of Hausdorff dimension p is uniformly distributed on the manifold for s ≥ p, also holds here.The case s = 2 deserves special attention.It is plausible that for some s we might have analytical solutions, for instance s = 2 (besides the usual s = 1 case).
After some considerable algebra, we arrive at the result ) r(r) where P s−2 2 (⋅) is the Legendre function of index s−2 2 .For s ≥ 2, we obtain that the only compatible result for (87) is that (r) has to remain constant.Mathematically, the integration over r is split between (0, r ′ ) and (r ′ , R).The ensuing result requires the use of the Cauchy principal value, which entails a uniform charge density.Now, the solution of (87) may be obtained by inverting the integral equation.This is an arduous task.On the other hand, one can show that a certain form for (r) reduces the right-hand side of (87) to a constant.The equilibrium density, once normalized, is given by with s ∈ (0, 2).The number of particles up for r ≤ R is given by . There is a brilliant analogy in contact problems for an elastically inhomogeneous half-space with a variable elastic modulus E(z) that goes as z k , for −1 < k < 1, which provides the corresponding Boussinesq problem of the fundamental problem of the kind w(r) ∝ 1∕r 1+k , or 1∕r s with s ∈ (0, 2).Thus, a varying exponent in the electrostatic case is to be found as a variable (vertically) elastic modulus for the inhomogeneous half-space.What is done there avoids the kind of integration in (87), as it is worth to be shown in detail.

Limiting Equilibrium Density for Charges Confined Inside the Sphere
Let us suppose that N charges are interacting via the modified Coulomb potential 1∕r s .The analytic expression for the equilibrium distribution when N → ∞ is well known; it appears, for example, in ref. [79] for a general dimension p.According to Hardin and Saff, [78] the distribution shall be uniform for s ≥ 3, which is supported in ref. [77] for p = 3.The corresponding variational calculation in dimension p leads to The limiting density then takes the form ) 2 ] p−s 2 (90) for s ∈ (p − 2, p), while (r) will be constant for s ≥ p. Equation (90) matches the previous result for s = 2 from ref. [77].Note that in the range −2 < s ≤ p − 2, the equilibrium is achieved for the uniform measure on the boundary of the ball-that is, on the sphere of radius R.This is in line with the property of classical electrical charge (corresponding to s = 1, p = 3) to concentrate on the surface of a conductor.

Summary and Conclusions
Potential theory is able to mathematically unify different parts of physics that apparently bear no connection, such as equilibrium electrostatics, heat transfer, or elasticity theory (contact mechanics).In the present contribution we aim to elucidate a more physical twist to the underlying formalism of potential theory for scientists interested in several research fields of mesoscopic systems dealing with charges in equilibrium.Although the connection between contact mechanics and electrostatics is known to the specialist, it is not so to the broader audience.Systems of confined charges treated quantum mechanically are often compared with their classical counterpart in the thermodynamic limit.Most of the existing results corresponding to several confining potentials have their equilibrium charge distributions scattered in the literature.Therefore, this work intends to fill this information gap.
In the process of comparing two apparently disjoint theories, we have been able to address the cases where several confin-ing potentials induce equilibrium analytic surface densities.Furthermore, we have introduced the tools-either numerical or analytical-usually employed in contact mechanics and in the framework of equilibrium electrostatics.Specifically, we have retrieved a series of equilibrium charge distributions for commonly used confinements.They range from linear radial ones to polar confinements, considering the cases which can be of more interest to the interested audience.
We have not limited ourselves to detail the correspondence between contact mechanics and equilibrium electrostatics for the sake of employing known results from the former applied to the latter.On the contrary, and motivated by the variational approach to finding the equilibrium charge distribution axisymmetric systems (particularly those with hard-wall confinements), we have studied systems not previously considered.The axial symmetry covers a considerable range of problems, and we have provided the analytic and numerical tools to analyze them.The charged cylinder without caps has been solved analytically and numerically, and the analytic solution for the charged cone has been tantalizing.
The case of the disk with two species has been proposed, but its solution is far from trivial.Also, the important case of the circular capacitor has been solved via the corresponding integral equation, and some discrete results have been shown.72] The self-energy for compact charged objects has also been provided, for the sake of generality.Also, variations of the Coulomb interaction in the form of 1∕r s have been presented with details.This particular interaction-the Riesz potential, is well known as an integral kernel for circular systems in many dimensions in potential theory.Surprisingly, the concomitant analytic results are mostly unknown to the physics community.Summing up, by extending the formalism of confined charges to axisymmetric 3D surfaces, specific examples have been worked in detail, providing new results previously unknown in the literature.

Figure 1 .
Figure 1.Graphical description of the Boussinesq problem.As defined in Podio-Guidugli and Favata,[37] "the Boussinesq problem (Joseph Valentin B., 1842-1929) consists in finding the elastic state in a linearly elastic isotropic half-space, subject to a concentrated load applied in a point of its boundary plane and perpendicular to it."An indenter on an half-infinite elastic media exerts stresses normally (pressure − zz and perpendicularly  zr ) to the point of action.The mathematical dependence on r of the action at a distance r perpendicular to the indentation is what makes the problem relevant mathematically.See text for details.

Figure 2 .
Figure2.Graphical analogy between axisymmetric contact mechanics and equilibrium electrostatics.a) The 3D profile of the rigid indenter f (r) that will punch an semi-infinite elastic media, which plays the role of the confining potential V(r) (up to a constant) for the continuum of charges.b) The action of a cylindrical flat punch on an elastic media induced by a force on top of the indenter, which is the typical profile studied in contact mechanics.c) The equivalent role played by a hard-wall confinement V(r).Visually, the analogy is perfect.d) The simplest non-trivial case after the flat punt is shown by the conical indenter, that corresponds to a linear confinement potential V(r). g

Figure 3 .
Figure 3. Distribution of charges in equilibrium for the triangle (N = 120), square (N = 1000), rectangle (N = 2000), and hexagon (N = 1200).In each case, a numerical optimization of the total interaction energy was obtained by an interior point optimization algorithm[62] with inequality constraints.When the number of charges N ≫ 1, one can observe a flat and compact distribution in the innermost regions, which implies an almost uniform charge density away from the boundaries.Close to the boundary, on the other hand, accumulation of charges indicates a spike in density.See text for details.

Figure 4 .
Figure 4. Equilibrium positions for N = 600 self-interacting charges on the lateral surface of the cylinder.The dimensions of the confining cylinder correspond to L = R = 1.See text for details.

Figure 5 .
Figure 5. Cone integral equation numerical solution for several ratios H∕R.The inset displays the equilibrium density for the hemisphere.See text for details.

Figure 7 .
Figure 7. Equilibrium on a capacitor versus a single disk.Left: Distribution of N = 1000 (on each plate) charges in equilibrium on the capacitor.Charges on different plates are superimposed, in blue and yellow.Right: Distribution of N = 1000 charges in equilibrium on a single disk, in red.Middle: Smoothed density histograms of the two types of distributions, capacitor (blue) and single disk (red).See text for details.

Figure 8 .
Figure 8. Capacitor integral equation numerical solution versus the distance between plates d∕R (in units of R) and the radius r∕R.See text for details.