For the calculation of the electrostatic potential U of a charge density one 
can use the equation (in atomic units)

Inserting the Laplace expansion

with the abbreviations , , and  
and the complex spherical harmonics in the phase convention of Condon and Shortley
 [126] as defined in [127], p. 3, Eq. (1.2-1), and interchanging the order of integration and 
summation yields

This can be interpreted as orthogonal expansion in the complete orthonormal system of 
spherical harmonics.

If the charge density vanishes outside of a sphere of radius a with center at the origin and 
if r>a holds then r>r' 
is satisfied for all for which holds, and the relation

follows.

Introduction of the multipole moments

allows to represent the potential by an infinite series of the form

that is usually called the multipole expansion. In particular, for a charge distribution that is
invariant under rotations with axis , the multipole moments have the form

and the addition theorem

of the spherical harmonics yields an expansion in Legendre polynomials

This form of the electrostatic potential is also implied by the fact that the Laplace 
equation  holds outside of a charge distribution. The only solutions of 
this equation for in spherical coordinates are linear combinations of irregular 
solid harmonics and thus, taking rotational symmetry 
around the direction into account, only linear combinations of  
occur in the multipole expansion of U where is the angle between the position vector and 
the direction . 
If the charge distribution does not have compact support in , the above multipole expansions 
only hold approximately for large . In this case, the inequality r>r' does not hold for 
all for which . 
The difference between the exact potential U and the multipole expansion in Eq. (6) is given by

Because of

the difference can also be represented as

If decays sufficiently rapidly for large arguments, the contributions with angular 
momentum quantum number to the difference will go to zero for large r  
(unless they vanish anyway by symmetry). In the case of the example
  with , only the term with survives the 
angular integration and one obtains 

Hence, the difference vanishes for large r exponentially in this simple 
model example.
Note that the multipole expansion for is a solution of the Laplace equation . 
But the exact electrostatic potential U satisfies the Poisson equation (atomic units)

This implies that also the difference is a solution of this Poisson equation. 
A further consequence  is that the multipole expansion can only be a good approximation 
to U where the charge density is small, i.e., for large distances from the origin.

An advantage of the multipole expansion is that the moments or , respectively, to a 
given charge distribution  can be computed once and for all, and then, the multipole 
approximation of can be computed very easily for very many arguments . 
The exact expansion (3) of on the other hand is more demanding computationally. 
It requires the calculation of the integrals 

that depend also on the distance r to the expansions centre. If the charge density is rotationally 
symmetric around the direction , then in analogy to Eq. (7) we have

and hence, the expansion (3) simplifies to

This corresponds to the cylindrical symmetry of the problem since only the coordinates r  
and with enter. Thus, one obtains an expansion in Legendre 
polynomials that may or may not converge rapidly.

Information provided by: http://www.chemie.uni-regensburg.de