A point charge is a point particle[N1A point particle defines a particle which lacks spatial
extension.] which has an electric charge. There are two types
of electric charges: positive and negative.
The elementary charge, a fundamental physical constant usually denoted as
e, is the electric charge carried by a single proton, whereas a single
electron has charge −e. In the International System of Units (SI)
the elementary charge has a measured value of approximately
1.602 176 6208(98)×10^{−19} C.[1(a) Mohr, P. J.; Newell, D. B.; Taylor, B. N.
Reviews of Modern Physics 2016, 88, 035009, 73 pages.
(b) Mohr, P. J.; Newell, D. B.; Taylor, B. N.
J. Phys.
Chem. Ref. Data 2016, 45, 043102, 74 pages.]
In the natural system of Hartree atomic units (usually abbreviated
“au” or “a.u.”)[N2Definition of Hartree atomic units.] the
elementary charge functions as the unit of electric charge, that is e
is equal to 1.
All observable charges occur in integral amounts of the
fundamental unit of charge e, that is, charge is quantized.
Any observable charge occurring in nature can be written
Q = ±Ne,
where N is an integer.[2Tipler,
Paul A.; Mosca, Gene
"Physics for Scientists and Engineers",
 5th ed.
W. H. Freeman and Company, San Francisco, 2004; Vol. I,
p. 653.]
The amount of electric charge per unit volume is
given by the charge density, ρ(r), whose SI units are
C⋅m^{−3}. The atomic unit of charge density
(e/a_{0}^{3}, where a_{0} = 0.529 177 210 67(12) × 10^{−10}
m is the radius of the first Bohr orbit for hydrogen) is equal to 1.081 202 3770(67) × 10^{12} C ⋅
m^{−3}.[1(a) Mohr, P. J.; Newell, D. B.; Taylor, B. N.
Reviews of Modern Physics 2016, 88, 035009, 73 pages.
(b) Mohr, P. J.; Newell, D. B.; Taylor, B. N.
J. Phys.
Chem. Ref. Data 2016, 45, 043102, 74 pages.]
For an ensemble of N point charges q_{i} located at points R_{i} (i = 1, 2, …, N), the total charge Q is the sum of the N individual point charges:
Q = N ∑ i=1 q_{i}  (1.1) 
and the total charge density is the sum of the N individual charge densities
ρ(r) = N ∑ i=1 q_{i} δ(r − R_{i})  (1.2) 
which are expressed in terms of the Dirac “delta” functions,
δ(r −
R_{i}).[3Dirac, P. “The Principles of Quantum Mechanics”
Oxford at the Clarendon Press: 1958.,
4Messiah, A. “Quantum Mechanics”,
Volume I
NorthHolland Publishing Company, Amsterdam: 1970,
pp. 179184 and 468471.,
N3ThreeDimensional Dirac Delta
Function.]
The function δ(r −
R_{i}) is equal to zero everywhere but at
the single point
r = R_{i}, where it
is singular in such a manner that
∫ δ(r − R_{i}) d^{3}r = 1  (1.3) 
Integrating eq. (1.2) over any volume that contains all the R_{i} yields, because of eq. (1.3), the total charge of the distribution, i.e. Eq (1.1)
Q = ∫ ρ(r) d^{3}r = N ∑ i=1 q_{i} ∞ ∫ −∞ δ(r − R_{i}) d^{3}r = N ∑ i=1 q_{i} 
Strictly speaking, it would not be possible to define a density function for a discrete distribution of charges, since the density is intrinsically a continuous function. The use of the Dirac “delta” function allows defining a generalized density function that substantially unifies the treatment of discrete and continuous distributions. This means, for instance, that formulas given for a continuous distribution of charges can be used as well for discrete distributions.
The sum of the vectors q_{i} R_{i} defines the dipole moment vector of the distribution
p = N ∑ i=1 q_{i} R_{i}  (1.7) 
In particular, an ensemble of two pointlike particles carrying charges of the same strength q and opposite sign represents an electric dipole. An electric dipole can be characterized by its electric dipole moment, a vector quantity that points from the negative charge (q_{1} = −q at location r_{1} = r_{−}) towards the positive charge (q_{2} = +q at location r_{2} = r_{+}), so that it is defined by the difference
p = q (r_{+} − r_{−})  (1.8) 
The magnitude of the electric dipole moment is given by the value of the charge q times the distance of the opposite charges. It is a measure of the polarity of the charge distribution, that is of the degree of local charge separation. In other words, a dipole exists when the centers of positive and negative charge distribution do not coincide. Though the position vectors r_{+} and r_{−} are defined with respect to some reference system of coordinates, it is clear from eq. (1.8) that the electric dipole moment of a pair of opposite charges does not depend on the choice of origin. This fact holds also in the general case of eq. (1.7), provided the overall charge of the system is zero. Indeed, if all the position vectors, r_{i}, are moved to a new origin r_{0} (i.e. r_{i} → r_{i} − r_{0}), then eq. (1.7) becomes:
p(r_{0}) = N ∑ i=1 q_{i} (r_{i} − r_{0}) = N ∑ i=1 q_{i} r_{i} − Q r_{0}  (1.9) 
which is independent of r_{0} only if Q = 0.[N?4Conventional choice of the origin for charged systems.]
The SI units for electric dipole moment are coulombmeter (C⋅m), however the most commonly used unit is the debye (D), which is a CGS unit. A debye is 10^{−18} esu⋅cm and about 3.336 × 10^{−30} C⋅m. An atomic unit of electric dipole moment is e⋅a_{0} and is about 8.478 × 10^{−29} C⋅m.
A general distribution of electric pointcharges may be characterized, in addition to its net charge and its dipole moment, by its electric quadrupole moment and, eventually, by higher order moments. The electric quadrupole moment is a symmetric ranktwo tensor (i.e. a 3×3 matrix with 9 components, but because of symmetry only 6 of these are independent).
Θ =
N
∑
i=1
q_{i} r_{i}
r_{i}^{†} =
N
∑
i=1
q_{i}
 (1.10) 
As with any multipole moment, if a lowerorder moment (monopole or dipole in this case) is nonzero, then the value of the quadrupole moment depends on the choice of the coordinate origin. If all the position vectors, r_{i} in eq. (1.10), are moved to a new origin r_{0} (i.e. r_{i</sub> → ri − r0), eq. (1.10) becomes:}
Θ(r_{0}) =
N
∑
i=1
q_{i} (r_{i} − r_{0})
(r_{i} −
r_{0})^{†} =
N ∑ i=1 q_{i} r_{i} r_{i}^{†} − p r_{0}^{†} − r_{0} p^{†} + Q r_{0} r_{0}^{†}  (1.11) 
which states that Θ is independent of r_{0} (i.e. does not depend on the choice of origin) only if p = 0 and Q = 0. For example, the quadrupole moment of an ensemble of four samestrength pointcharges, arranged in a square with alternating signs, is coordinate independent.
The SI units for electric quadrupole moment are C⋅m^{2}, however the most commonly used unit is the buckingham (B), which is a CGS units. A buckingham is 1 × 10^{−26} statC⋅cm^{2} and is equivalent to 1 D⋅Å. An atomic unit of electric quadrupole moment is e⋅a0^{2} and is about 4.486 551 484(28) × 10^{−40} C⋅m^{2}.
“Often it is convenient to ignore the fact that charges come in packages like electrons and protons, and think of them as being spread out in a continuous smear − or in a 'distribution,' as it is called. This is O.K. so long as we are not interested in what is happening on too small a scale. We describe a charge distribution by the 'charge density,' ρ(r) ≡ ρ(x,y,z). If the amount of charge in a small volume ΔV_{1} located at the point r_{1} is Δq_{1}, then ρ is defined by”[L1“Coulomb’s law; superposition”
The Feynman Lectures on Physics, Online Edition,
Volume II, Chapter 4, Section 4.2]
Δq_{1} = ρ(r_{1}) ΔV_{1}  (2.1) 
As always, the volume integral of the charge density ρ(r) over a region of space V in ℝ^{3} is the charge contained in that region.
q = ∫ V ρ(r) d^{3}r = ∫ V ρ(r) dV = ∭ V ρ(x,y,z) dx dy dz  (2.2) 
In the Euclidean space, the components of the position vector r are the cartesian coordinates x, y and z, and the volume element dV is given by the product of the differentials of the cartesian coordinates, dV = dx dy dz.
Matter consists of atoms that are electrically neutral.
Each atom has a tiny but massive nucleus that contains protons and neutrons.
Protons are positively charged, whereas neutrons are uncharged.
The number of protons in the nucleus is the atomic number Z of the
element. Surrounding the nucleus is an equal number of negatively charged
electrons, leaving the atom with zero net charge.
The electron is about 2000 times less massive than the proton, yet the
charges of these two particles are exactly equal in magnitude.[2Tipler, Paul A.; Mosca, Gene
"Physics for Scientists
and Engineers",  5th ed.
W. H. Freeman and Company, San Francisco, 2004;
Vol. I, p. 653.]
The Heisenberg's uncertainty principle states that the
more precisely the position of an electron is determined, the less precisely
its momentum can be known, and vice versa.[5Heisenberg, W. Zeitschrift für Physik
1927, 43, 172–198., 6Sen, D. Current Science 2014, 107,
203218.]
Therefore, we can not say exactly where the electron is, but quantum
mechanics provides a recipe for calculating the probability
ρ(x,y,z)
dx dy dz that it is in an element of volume
dV = dx dy
dz at a particular point in the coordinate space x,
y, z, and this quantity is expressed in terms of a function
ρ(x,y,z)
called “electron probability density” or simply “electron
density”.
This does not imply that the electron is “diluted” throughout
the region of space in which the electronic density is nonzero: in a certain
volume dV there may be a very small probability of finding the electron; but if you find it there, you find it all.
However, in many cases (for example, dispersion of Xrays, forces on atoms),
a system consisting of a large number of electrons behaves exactly as if the
electrons were distributed nonuniformly over a region of space, so that each
point in space is characterized by a certain fraction of charge or partial
charge, which is identical to the electron probability density.
In fact, for large numbers, the frequency tends to probability.
Therefore, a map showing how the electron probability density varies from
point to point in a molecule, i.e. a map of the cumulative distribution of
the charge of all the electrons of the molecule, is a model of the
electronic structure of the molecule.[7Shusterman, G. P.; Shusterman, A. J.
J. Chem. Ed.
1997, 74, 771776.]
Two or more atoms held together by chemical bonds give an electrically neutral entity called molecule.[8IUPAC. Compendium of Chemical Terminology, 2nd ed. (the "Gold Book").] Since atomic nuclei are much heavier than electrons (about 2000 times), they move more slowly. Hence, to a good approximation (the BornOppenheimer approximation[9Born, M.; Oppenheimer, R. Annalen der Physik 1927, 389, 457484.]), we can consider that the electrons are moving in a field of fixed nuclei.
The molecular charge density ρ_{molecule}(r,{R_{i}}) is the sum of the nuclear density ρ_{nuclei}(r,{R_{i}}) and the electron density ρ_{electrons}(r,{R_{i}}):
ρ_{molecule}(r,{R_{i}}) = ρ_{nuclei}(r,{R_{i}}) − ρ_{electrons}(r,{R_{i}})  (3.1) 
where the variable r is the position vector and the nuclear coordinates {R_{i}} are fixed parameters. The − sign accounts for the negative charge of the electrons.
The nuclei in a molecule are a discrete distribution of positive charges {Z_{i}} at fixed positions {R_{i}} so that their density has the following expression:
ρ_{nuclei}(r,{R_{i}}) = N ∑ i=1 Z_{i} δ(r − R_{i})  (3.2) 
according to eq. (1.2). The charge of each nucleus is given by its atomic number Z_{i}. The expression of the electron density, a continous probability distribution, depends on the quantummechanical recipe used for its calculation or by the functional model used to analyze experimental xray diffraction data. Eventually, it can be known as a set of numerical values defined over a grid of data points {r_{j}}.
Any molecular probability charge distribution involves both a continuous and a discrete part (eq. 3.1). As noted previously, the use of the generalized probability density function (eq. 3.2) substantially unifies the treatment of continuous and discrete distributions.
4.1.  Joint and Marginal Probability Density Function. In the previous section, the electron probability density has been introduced as the scalarvalued function ρ(x,y,z) of the three Cartesian coordinates x, y, z in the Euclidean space ℝ^{3} such that, multiplied by an infinitesimal element of volume dV = dx dy dz at point (x,y,z), gives the probability of finding the electron at that point. This definition can be extended to any charge distribution.
More correctly, the function ρ(x,y,z) should be called a joint probability density function, where the word “joint” stresses the fact that this probability is linked to the occurrence of a triple of values, one for each of the three variables x, y, z. In common practice, it is simply called probability density function (PDF), or even more shortly, probability function or density function.
The PDF has the following properties:
A number of marginal PDF's can be obtained from the joint PDF ρ(x,y,z) by means of the following definitions:
ρ(x) = ∞ ∫ −∞ ∞ ∫ −∞ ρ(x,y,z) dy dz  (4.1.1) 
ρ(y) = ∞ ∫ −∞ ∞ ∫ −∞ ρ(x,y,z) dx dz  (4.1.2) 
ρ(z) = ∞ ∫ −∞ ∞ ∫ −∞ ρ(x,y,z) dx dy  (4.1.3) 
ρ(x,y) = ∞ ∫ −∞ ρ(x,y,z) dz  (4.1.4) 
ρ(x,z) = ∞ ∫ −∞ ρ(x,y,z) dy  (4.1.5) 
ρ(y,z) = ∞ ∫ −∞ ρ(x,y,z) dx  (4.1.6) 
Therefore, the PDF associated with each variable x,
y, and z alone (or with any pair of variables xy,
xz, and yz) is called marginal density function,
and can be deduced from the joint probability density associated with the
variables x,y,z by integrating over all values of the remaining two
(or one) variables.
The radial density function is the most important example of
marginal PDF that we will find in this manual.
It proves useful to analyse the PDF with the aid of related functions which point out features and properties of the PDF itself. A few of them will be considered in the following.
4.2.  Cumulative Distribution Function (CDF). The (joint) cumulative distribution function D(x,y,z) is defined as the probability to find some charge within the semiinfinite parallelepiped with upper limits defined by x, y, z. It can be expressed as the integral of its (joint) probability density function ρ(x,y,z) as follows:
D(x,y,z) = x ∫ −∞ y ∫ −∞ z ∫ −∞ ρ(x',y',z') dx' dy' dz'  (4.2.1) 
The (joint) cumulative distribution function D(x,y,z) has the following properties:
Of course, a CDF exists for all the marginal PDFs of ρ(x,y,z), i.e.:
D(x,y) = x ∫ −∞ dx' y ∫ −∞ dy' ∞ ∫ −∞ ρ(x',y',z') dz' = x ∫ −∞ y ∫ −∞ ρ(x',y') dx' dy'  (4.2.2) 
D(x,z) = x ∫ −∞ dx' z ∫ −∞ dz' ∞ ∫ −∞ ρ(x',y',z') dy' = x ∫ −∞ z ∫ −∞ ρ(x',z') dx' dz'  (4.2.3) 
D(y,z) = y ∫ −∞ dy' z ∫ −∞ dz' ∞ ∫ −∞ ρ(x',y',z') dx' = y ∫ −∞ z ∫ −∞ ρ(y',z') dy' dz'  (4.2.4) 
D(x) = x ∫ −∞ dx' ∞ ∫ −∞ ∞ ∫ −∞ ρ(x',y',z') dy' dz' = x ∫ −∞ ρ(x') dx'  (4.2.5) 
D(y) = y ∫ −∞ dy' ∞ ∫ −∞ ∞ ∫ −∞ ρ(x',y',z') dx' dz' = y ∫ −∞ ρ(y') dy'  (4.2.6) 
D(z) = z ∫ −∞ dz' ∞ ∫ −∞ ∞ ∫ −∞ ρ(x',y',z') dx' dy' = z ∫ −∞ ρ(z') dz'  (4.2.7) 
In turn, the PDFs can be obtained as the derivative of their CDF:
ρ(x) = dD(x) dx  (4.2.8) 
ρ(y) = dD(y) dy  (4.2.9) 
ρ(z) = dD(z) dz  (4.2.10) 
ρ(x,y) = ∂^{2}D(x,y) ∂x ∂y  (4.2.11) 
ρ(x,z) = ∂^{2}D(x,z) ∂x ∂z  (4.2.12) 
ρ(y,z) = ∂^{2}D(y,z) ∂y ∂z  (4.2.13) 
ρ(x,y,z) = ∂^{3}D(x,y,z) ∂x ∂y ∂z  (4.2.14) 
In their 1968 study of the four lowest
^{1}Σ^{+} states of the LiH molecule, Brown and Shull
used the cumulative distribution function D(z) of Eq. (4.2.7)
and the marginal probability density function ρ(z) of
Eq. (4.2.10) to analyse the calculated wave functions.[10Brown, Richard E.; Shull, Harrison
Int. J.
Quantum Chem. 1968, 2, 663685.]
These two functions were designated respectively electron count and
planar density. Since the Li and H nuclei are situated on the
zaxis with the axis origin at the molecular midpoint, the electron
count function D(z) simply denotes the number of electrons on
the negative side of the xyplane intersecting the zaxis at
the point z, and the planar density
ρ(z) = D'(z)
is the density of electrons in this plane.[10Brown, Richard E.; Shull, Harrison
Int. J.
Quantum Chem. 1968, 2, 663685.]
In 1970, Politzer and Harris used the electron count function D(z)
to propose a definition of the charge on an atom in linear molecules.[11Politzer, Peter; Harris, Roger R.
J.
Am. Chem. Soc. 1970, 92, 64516454.]
In 1979 Streitwieser and coworkers suggested to use the
bivariate marginal density ρ(x,z) of Eq. (4.5) or
(4.2.12), called the electron projection function, for aiding
visualization of the electron density.
Projecting the electron density function ρ(x,y,z) into
the xzplane corresponds to summing
ρ(x,y_{i},z) over all
y_{i} = constant
parallel planes. Hence, there is no arbitrariness about the choice of one of
these planes in which to display the PDF.
Because it compresses the total electron density onto a twodimensional grid,
it is ideal for graphical representation and simplifies qualitative electron
density analyses.[12(a)
Streitwieser, A. Jr.; Collins, J. B.; McKelvey, J. M.; Grier, D.; Sender,
J.; Toczko, A. G.
Proc. Nati. Acad. Sci. USA 1979, 76,
24992502.
(b) Collins, J. B.; Streitwieser, A. Jr.
J. Comput. Chem. 1980, 1, 8187.
(c)
Bachrach, S. M.; Streitwieser, A. J. Comput. Chem. 1989,
10, 514519.
(d) Agrafiotis, D. K.; Tansy, B.;
Streitwieser, A. J. Comput. Chem. 1990, 11,
11011110.]
The electron projection function was also proposed as a generalization for
planar or nearly planar systems of the partitioning of space introduced by
Politzer[11Politzer, P.; Harris, R. R.
J. Am. Chem. Soc. 1970, 92, 64516454.] for
linear molecules. However the projection method gives integrated charges that
have no fundamental significance and which are at best only approximations to
charges that correspond to zeroflux surfaces.[13Glaser, R. J. Comput. Chem. 1989, 10,
118135.]
4.3.  Derivatives. The density function ρ(x,y,z) = ρ(r) is a continuous function which operates on a vector variable r = (x, y, z)^{†} ∈ ℝ^{3} and returns a scalar, i.e. ρ(r): ℝ^{3} → ℝ. Since ρ(r) assigns a scalar value to any point r in space, it defines a scalar field. The function can be differentiated to any order with respect to each of the variables x, y, and z, the other variables being understood to be held constant during the differentiation. Such derivatives are known as partial derivatives. The firstorder partial derivatives are
∂ ρ ∂ x = ρ_{x}, ∂ ρ ∂ y = ρ_{y}, ∂ ρ ∂ z = ρ_{z}  (4.3.1) 
As the variables x, y, z are the components of the (column) vector r, similarly ∂ ∂ x , ∂ ∂ y , ∂ ∂ z are the components of the gradient vector operator, “grad”, more commonly indicated by the symbol ∇, denominated “nabla” or “del”.[N4The vector differential operator ∇.] So, the three scalar expressions of Eq (4.3.1) are equivalent to the vector expression
grad ρ = ∇ ρ =
∂
ρ(r)
∂
r
=
 (4.3.2) 
The gradient of ρ(r) is a mapping which operates on a continuous scalar function of a vector variable and returns a continuous vectorvalued function, i.e. ∇ ρ(r): ℝ^{3} → ℝ^{3}. It is a physical quantity that has both magnitude and direction. Since it assigns a vector to any point r in space, it defines a vector field. Its components are organized in a column, so that it is referred to as a column vector or column matrix, i.e. a matrix which has only one column. The transpose of a column matrix is another matrix which has only one row and is called a row matrix or row vector.[N5Definition of column matrix and column vector]
By differentiating the three firstorder partial derivatives of Eq. (4.3.1) with respect to x, y, and z, three secondorder partial derivatives are obtained
∂^{2} ρ ∂ x^{2} = ∂ ρ ∂ x ( ∂ ρ ∂ x ) = ρ_{xx}, ∂^{2} ρ ∂ y^{2} = ∂ ρ ∂ y ( ∂ ρ ∂ y ) = ρ_{yy}, ∂^{2} ρ ∂ z^{2} = ∂ ρ ∂ z ( ∂ ρ ∂ z ) = ρ_{zz}  (4.3.3) 
together with six secondorder mixed derivatives
ρ_{xy} = ∂^{2} ρ ∂ x ∂ y = ∂ ρ ∂ x ( ∂ ρ ∂ y ) = ∂ ρ ∂ y ( ∂ ρ ∂ x ) = ∂^{2} ρ ∂ y ∂ x = ρ_{yx},  
ρ_{xz} = ∂^{2} ρ ∂ x ∂ z = ∂ ρ ∂ x ( ∂ ρ ∂ z ) = ∂ ρ ∂ z ( ∂ ρ ∂ x ) = ∂^{2} ρ ∂ z ∂ x = ρ_{zx},  (4.3.4) 
ρ_{yz} = ∂^{2} ρ ∂ y ∂ z = ∂ ρ ∂ y ( ∂ ρ ∂ z ) = ∂ ρ ∂ z ( ∂ ρ ∂ y ) = ∂^{2} ρ ∂ z ∂ y = ρ_{zy} 
where it appears that the mixed derivatives (assuming that they are continuous) do not depend on the order of the differentiation. Higherorder partial and mixed derivatives are again defined recursively, e.g.
∂^{i+j+k} ρ ∂x^{i} ∂y^{j} ∂z^{k} = ρ^{(i,j,k) }  (4.3.5) 
The second order partial and mixed derivatives of Eqs (4.3.3) and
(4.3.4) are found in the result of applying the “nabla”
operator twice in sequence. How they are structured in the result
depends on how the “nabla” operator is applied.
Since ∇ can be viewed as a 3×1 matrix,[N5Definition of column matrix and column
vector] let's take the Kronecker product[14Zhang, H.; Ding, F.
Journal of Applied
Mathematics 2013, Article ID 296185, 8 pages.,
N6Kronecker product
(definition),
L4Wikipedia contributors, "Kronecker
product,"
Wikipedia, The Free Encyclopedia.,
L5Weisstein, E. W. "Kronecker
Product."
From MathWorld  A Wolfram Web Resource.]
of ∇ with its transpose ∇^{†}
∇⊗∇^{†}ρ =
 (4.3.6) 
If all second partial derivatives of ρ(r)
exist and are continuous over the domain of the function, then the
Hessian matrix H of ρ(r) is a
square symmetric 3×3 matrix, defined and arranged as in
Eq (4.3.6).[L6Wikipedia contributors,
"Hessian matrix,"
Wikipedia, The Free Encyclopedia.]
The Kronecker product
∇^{†}⊗∇ =
(∇⊗∇^{†})^{†}
can be calculated in a similar way, leading to the definition of the
Jacobian matrix:[L7Wikipedia
contributors, "Jacobian matrix and determinant,"
Wikipedia, The Free
Encyclopedia., L8Weisstein,
E. W. "Jacobian."
From MathWorld  A Wolfram Web
Resource.]
∇^{†}⊗∇ ρ
=
 (4.3.7) 
The Jacobian matrix is the matrix of all firstorder partial
derivatives of the vectorvalued function
∇ρ(r), defined and arranged
as in Eq. (4.3.7). (Note that some literature defines the Jacobian
as the transpose of the matrix given in Eq. (4.3.7).) It defines a
linear map ℝ^{3} → ℝ^{3}, which is
the generalization of the usual notion of derivative, and is called
the derivative or the gradient of ∇ρ
at r.[L9Granados, A.
"Taylor Series For MultiVariable Functions." 2015.
A ResearchGate Web Resource.]
The Jacobian matrix of the gradient of the scalar function ρ(r) is the transpose of the Hessian matrix, and since the latter is a symmetric square matrix, they are equal. In a sense, both the gradient and Jacobian are “first derivatives”: the former the first derivative of a scalar function of r, the latter the first derivative of a vector function of r. The Hessian matrix in a sense is the “second derivative” of a scalar function of r.
H[ρ(r)] = J[∇ρ(r)]^{†}  (4.3.8) 
The Kronecker product ∇⊗∇ supplies the same derivatives of the Jacobian matrix in Eq. (4.3.7), but organized in a single 9×1 column matrix:
∇⊗∇ ρ =
 (4.3.9) 
where the “vec” (vectorizaton) operator has been
introduced.[15Magnus, J. R.; Neudecker,
H.
The Annals of Statistics 1979, 7,
381394., 16Henderson, H. V.;
Searle, S. R.
Linear and Multilinear Algebra 1981,
9, 271288., L10Wikipedia contributors, "Vectorization_(mathematics),"
Wikipedia, The Free Encyclopedia.]
It is worth noting that the same result of Eq. (4.3.6) can be obtained by multiplying ∇ and ∇^{†} according to the usual matrix multiplication rule, provided that ∇ is represented as a 3×1 column matrix (which makes ∇^{†} a 1×3 row matrix):
∇∇^{†}ρ =
 (4.3.10) 
i.e. ∇⊗∇^{†} = ∇∇^{†}. At this point we have the rules to evaluate all the thirds and higher derivatives of ρ(r) at any point r = (x, y, z)^{†} ∈ ℝ^{3}. For example, the derivative or the gradient of the Jacobian matrix J(r) is the 3×9 matrix of the thirds derivatives of ρ(r), defined and arranged as follows:
∇^{†}⊗J =
 (4.3.11) 
Similarly, the fourth derivatives of ρ(r) are the elements of the 3×27 matrix ∇^{†}⊗∇^{†}⊗J = ∇^{⊗2}J = ∇^{⊗3}∇ ρ(r) i.e. the fourth derivatives of ρ are the second derivatives of the Jacobian matrix and the third derivatives of the gradient of ρ.
In general, we'll use the notation of Eq. (4.3.9) to compute the derivatives of the scalar field ρ(r), so that its nth derivatives are the components of the 3^{n}×1 column matrix of Eq. (4.3.12):
∇⊗∇⊗ … ∇⊗ ntimes ρ(r) = ∇^{⊗(n)} ρ(r) = ∇^{⊗(n−1)}⊗∇ ρ(r) = ∇^{⊗(n−2)}⊗ J(r)  (4.3.12) 
where sometimes the symbol ⊗ is omitted.
The dot product of the vector operator ∇ with the gradient vector of ρ(r) defines the Laplacian of ρ(r), a scalar field:
∇⋅∇ρ = ∇^{2}ρ =
Δρ =
 (4.3.13) 
It's worth noting that the Laplacian of ρ(r) is equal to the trace of the Hessian (Jacobian) matrix of ρ(r):
Δρ(r) = ∇^{2}ρ(r) = Tr H(r) = Tr J(r)  (4.3.14) 
The Laplacian can also be defined as the divergence of the gradient of ρ(r).
The directional derivative of the scalar function ρ(r) = ρ(x,y,z), differentiable at r, along a vector u = (u_{x}, u_{y}, u_{z})^{†} is the scalar function ∇_{u} ρ(r) defined by the dot product of u with the gradient of ρ(r):
∇_{u} ρ(r) = u ⋅ ∇ρ(r)  (4.3.15) 
If u is a unit vector, the directional derivative ∇_{u} ρ(r) is independent of the magnitude of u and depends only on its direction. Intuitively, the directional derivative of ρ at a point r represents the rate of change of ρ with respect to a change of its argument r in the direction of u.
he Laplacian Δf(p) of a function f at a point p, is (up to a factor) the rate at which the average value of f over spheres centered at p deviates from f(p) as the radius of the sphere grows.4.4.  Taylor series expansion. Suppose ρ(r): ℝ^{3} → ℝ is of class C^{ k+1} on an open convex set S [i.e. ρ(r) is a scalar function of a vector variable r = (x, y, z)^{†} ∈ ℝ^{3}, its partial derivatives exist up to the (k+1)th order and are continuous, it is defined over a set of points such that, given any two points A, B in that set, the line AB joining them lies entirely within that set]. If a ∈ S and r = a + h ∈ S, then the Taylor expansion for ρ(r) at point r about point a is
ρ(a + h) = k ∑ j=0 h^{⊗j} ⋅ ∇^{⊗j} ρ(a) j! + R_{a,k}(h)  (4.4.1) 
where h = r − a, and formulas for the
remainder R_{a,k}(h) can be obtained
from the Lagrange or integral formulas for remainders, but this is beyond
the aim of this manual.[17Stewart, James
"Calculus Early Transcendentals", 6th ed.
Thomson Brooks/Cole,
USA, 2008, Section 11.10, p. 734., L14Stewart, J.
"Formulas for the Remainder Term
in Taylor Series"]
To obtain Eq. (4.4.1), we start by supposing that ρ can be
represented by a Kronecker power series[17Stewart, James "Calculus Early Transcendentals",
6th ed.
Thomson Brooks/Cole, USA, 2008, Section 11.10,
p. 734.]
ρ(a + h) =
k
∑
j=0
c_{j} ⋅ h^{⊗j} = c_{0} + c_{1} ⋅ h + c_{2} ⋅ h^{⊗2} + c_{3} ⋅ h^{⊗3} + c_{4} ⋅ h^{⊗4} + … + c_{k} ⋅ h^{⊗k}  (4.4.2) 
Let’s try to determine what the coefficients c_{j} must be in terms of ρ. To begin, notice that if we put r = a in Eq. (4.4.2), then h = 0, and all terms after the first one are 0 and we get
ρ(a) = c_{0}  (4.4.3) 
Since ρ is of class C^{ k+1}, we can differentiate the series in Eq. (4.4.2) term by term:
∇ρ(a + h) = c_{1} + 2 c_{2} ⋅ h + 3 c_{3} ⋅ h^{⊗2} + 4 c_{4} ⋅ h^{⊗3} + … + k c_{k} ⋅ h^{⊗(k−1)}  (4.4.4) 
and substitution of h = 0 in Eq. (4.4.4) gives
∇ρ(a) = c_{1}  (4.4.5) 
Now we differentiate both sides of Eq. (4.4.4) and obtain
∇⊗∇ρ(a + h) = ∇^{⊗2}ρ(a + h) = 2 c_{2} + 2⋅3 c_{3} ⋅ h + 3⋅4 c_{4} ⋅ h^{⊗2} + … + (k−1)⋅k c_{k} ⋅ h^{⊗(k−2)}  (4.4.6) 
Again we put h = 0 in Eq. (4.4.6). The result is
∇^{⊗2}ρ(a) = 2 c_{2}  (4.4.7) 
Let’s apply the procedure one more time. Differentiation of the series in Eq. (4.4.6) gives
∇⊗∇⊗∇ρ(a + h) = ∇^{⊗3}ρ(a + h) = 2⋅3 c_{3} + 2⋅3⋅4 c_{4} ⋅ h + 3⋅4⋅5 c_{5} ⋅ h^{⊗2} + … + (k−2)⋅(k−1)⋅k c_{k} ⋅ h^{⊗(k−3)}  (4.4.8) 
and substitution of h = 0 in Eq. (4.4.8) gives
∇^{⊗3}ρ(a) = 2⋅3 c_{3} = 3! c_{3}  (4.4.9) 
By now you can see the pattern. If we continue to differentiate and substitute h = 0, we obtain
∇^{⊗(k)}ρ(a) = 2⋅3⋅4⋅…⋅k c_{k} = k! c_{k}  (4.4.10) 
Solving this equation for the kth coefficient c_{k}, we get
c_{k} = ∇^{⊗(k)}ρ(a) k!  (4.4.11) 
This formula remains valid even for k = 0 if we adopt the conventions that 0! = 1 and ∇^{⊗(0)}ρ = ρ. Thus we have proved that if ρ has a Kronecker power series representation (expansion) at point a, as shown in Eq. (4.4.2), then its coefficients are given by Eq. (4.4.11). Substituting expression (4.4.11) for c_{k} back into the series (4.4.2), we see that if ρ has a Kronecker power series expansion at a, then it must be of the form given by Eq. (4.4.1).
For the special case a = 0, the Taylor series of Eq. (4.4.1) becomes
ρ(r) = k ∑ j=0 r^{⊗j} ⋅ ∇^{⊗j} ρ(0) j! + R_{0,k}(r)  (4.4.12) 
This case arises frequently enough that it is given the special name Maclaurin series.
4.5.  Moments and Moment Generating Function (MGF). The nth rank moment of a discrete or continuous probability distribution, described by a probability density function ρ(r) is given by the volume integral
m^{(n)} = ∫ V ρ(r) r^{⊗n} d^{3}r = E[r^{⊗n}]  (4.5.1) 
where E is the expectation operator or mean. The discussion of these quantities (that can be called “unnormalized traced Cartesian” moments about zero and are usually referred to as raw or crude moments[Weisstein, Eric W. "Raw Moment." From MathWorldA Wolfram Web Resource. http://mathworld.wolfram.com/RawMoment.html] and unabridged moments[Applequist J. Chem Phys. 1984; 85:279–290]) is postponed to a specific page of the manual. Here, we simply note that the zeroth rank moment of the distribution is just the electric charge defined by Eq. (2.2).
Given a discrete or continuous probability distribution, described by a probability density function ρ(r), the moment generating function (MGF) is defined by
M_{r}(s) = E(e^{s†r})  (4.5.2) 
for s finite, where r = (x, y, z)^{†} ∈ ℝ^{3}, s = (s_{x}, s_{y}, s_{z})^{†} ∈ ℝ^{3}, and
E(e^{s†r}) ≡ ⟨e^{s†r}⟩ = ∫ ℝ^{3} e^{s†r} ρ(r) d^{3}r  (4.5.3) 
is the expectation value of the function e^{s†r} in the variable r. Note that M_{r}(−s) is the twosided (doubly infinite) Laplace transform of ρ(r), ℒ_{r}[ρ(r)](s).
As its name implies, the moment generating function can be used to compute the moments of ρ(r): the kth moment about 0 is the kth derivative of the moment generating function, evaluated at s = 0. Let k_{x}, k_{y}, k_{z} be nonnegative integers and k = k_{x} + k_{y} + k_{z}. Then
∂^{k} ∂s_{x}^{kx} ∂s_{y}^{ky} ∂s_{z}^{kz} M_{r}(s)  = ∂^{k} ∂s_{x}^{kx} ∂s_{y}^{ky} ∂s_{z}^{kz} E(e^{s†r}) 
= E ( ∂^{k} ∂s_{x}^{kx} ∂s_{y}^{ky} ∂s_{z}^{kz} e^{sxx + syy + szz} )  
= E(x^{kx} y^{ky} z^{kz} e^{s†r}) 
Setting s = 0 = (0, 0, 0)^{†}, we get
∂^{k} ∂s_{x}^{kx} ∂s_{y}^{ky} ∂s_{z}^{kz} M_{r}(s)_{s=0} = E(x^{kx} y^{ky} z^{kz})  (4.5.4) 
An intuitively attractive expression for the function M_{r}(s), which accounts for the name “moment generating function”, can be obtained by expanding the function e^{s†r} = e^{s⋅r} as a Taylor series:
M_{r}(s) = E[e^{s⋅r}] = E [ ∞ ∑ j=1 s^{⊗j} ⋅ r^{⊗j} j! ] = ∞ ∑ j=1 E[s^{⊗j} ⋅ r^{⊗j}] j! = ∞ ∑ j=1 s^{⊗j} j! ⋅ E[r^{⊗j}] = ∞ ∑ j=1 s^{⊗j} j! ⋅ m^{(j)}  (4.5.5) 
So, if M_{r}(s) exists in a neighbourhood of 0, then the kth derivative of M_{r}(s) evaluated at s = 0 gives the kth moment m^{(k)} of the distribution:
m^{(k)} = E[r^{⊗k}] = ∇^{⊗k} M_{r}(s)_{s=0}  (4.5.6) 
which is the equivalent of Eq. (4.5.4) in tensor notation.
4.6.  Characteristic function. The characteristic function of the random vector r = (x, y, z)^{†} ∈ ℝ^{3}, which has a probability density function ρ(r), is defined as the expectation value of e^{is⋅r}, where i = √−1 is the imaginary unit, and s = (s_{x}, s_{y}, s_{z})^{†} ∈ ℝ^{3} is the argument of the characteristic function φ_{r}: ℝ^{3} → ℂ,
φ_{r}(s)  = E[e^{is⋅r}] = ∫ ℝ^{3} e^{is⋅r} ρ(r) d^{3}r  (4.6.1)  
= E[cos(s⋅r) + i sin(s⋅r)]  
= E[cos(s⋅r)] + i E[sin(s⋅r)] 
where the last equality in the first line is the Fourier transform of the probability density function ρ(r) with sign reversal in the complex exponential. Observe that φ_{r}(s) exists for any s ∈ ℝ^{3} because the expectation values appearing in the last line are welldefined, as both the sine and the cosine are bounded (they take values in the interval [−1,1]). So e^{is⋅r} lies on the unit circle in the complex plane, so that e^{is⋅r} = 1 and giving that φ_{r}(s) ≤ 1 for all s ∈ ℝ^{3}. Instead, the finiteness of the moment generating function of Eq. (4.5.2) is questionable since the exponential function is unbounded, so that e^{s⋅r} may be very large indeed.
The characteristic function and the probability density function are each a Fourier transform of the other. So, if the characteristic function is known, the inversion theorem allows to calculate the probability density function:
ρ(r) = 1 2 π ∫ ℝ^{3} e^{−is⋅r} φ_{r}(s) d^{3}s  (4.6.2) 
By expanding the function e^{is⋅r} as a Taylor series, the characteristic function can be written in the form
φ_{r}(s) = E[e^{is⋅r}] = E [ ∞ ∑ j=1 i^{ j} s^{⊗j} ⋅ r^{⊗j} j! ] = ∞ ∑ j=1 i^{ j} E[s^{⊗j} ⋅ r^{⊗j}] j! = ∞ ∑ j=1 i^{ j} s^{⊗j} j! ⋅ E[r^{⊗j}] = ∞ ∑ j=1 i^{ j} s^{⊗j} j! ⋅ m^{(j)}  (4.6.3) 
The characteristic function can be used to find moments of the random vector r. Provided that the kth moment exists, the characteristic function can be differentiated k times, yielding
∇^{⊗k} φ_{r}(s) _{s=0} = i^{k} E[r^{⊗k}] = i^{k} m^{(k)}  (4.6.4) 
4.7.  Transformation of the Random Vector Variable. Suppose we are given a random vector variable r ∈ ℝ^{3} with probability density ρ(r). If u ∈ ℝ^{3} is a random vector variable, we define a vectorvalued function f: ℝ^{3} → ℝ^{3} which takes as input the vector u and produces as output the vector r. In other words, if r_{1} = x, r_{2} = y, r_{3} = z are Cartesian coordinates, and u_{1}, u_{2}, u_{3} are some other coordinates, then the functions f_{1}, f_{2}, and f_{3} establish a onetoone correspondence between the coordinate systems: they have to be continuous, have continuous partial derivatives and single valued inverse, i.e.
r = f(u) or
 (4.7.1) 
Eq. (4.7.1) defines the socalled usubstitution.
Differentiating the vector r = (x, y, z)^{†} = (f_{1}, f_{1}, f_{1})^{†} = f yields the vector dr = (dx, dy, dz)^{†} = (df_{1}, df_{2}, df_{3})^{†} = df, which can be written more explicitly and with matrix notation as
 (4.7.2) 
or more concisely
dr = J_{f}(u) du  (4.7.3) 
The 3×3 matrix in the last equality of Eq. (4.7.2) is the Jacobian matrix of the function f at a point u where f is differentiable. It is usually denoted by
J_{f}(u) = Df =
∂ (f_{1},
f_{2}, f_{3})
∂
(u_{1}, u_{2},
u_{3})
=
∂ f
∂
u
= ∇^{†}⊗f =

∂ f
∂
u_{1}
∂ f
∂
u_{2}
∂ f
∂
u_{3}

=
 (4.7.4) 
or, componentwise:
J_{ij} =
∂ f_{i}
∂
u_{j}.
Notations in the second, third, and fourth members of Eq. (4.7.4)
point out that the Jacobian matrix can be seen as the generalization of
the usual notion of derivative, and then it is called the derivative or
the differential of f at u.[L7Wikipedia contributors, "Jacobian matrix and
determinant,"
Wikipedia, The Free Encyclopedia.,
L8Weisstein, E. W. "Jacobian."
From MathWorld  A Wolfram Web Resource.]
In particular, the last three members of Eq. (4.7.4) show that the
Jacobian matrix can also be thought as the gradient of f.
The vector
∂ f
∂
u_{j}
is the jth column of the Jacobian matrix (it is tangent
to the surface f at point u), whereas
the transpose of vector
∇ f_{i} = grad(f_{i})
is the ith row of the Jacobian matrix.
As a mnemonic, it is worth noting that across a row of the Jacobian
matrix the numerators are the same and down a column the denominators
are the same.
Leftmultiplying vectors in both members of Eq. (4.7.3) by their transpose yields the following equalities:
dr^{†} dr = du^{†} J_{f}(u)^{†} J_{f}(u) du = du^{†} g du  (4.7.5) 
where the symmetric matrix g = J_{f}(u)^{†} J_{f}(u) is called the fundamental (or metric) tensor of the Euclidean space in curvilinear coordinates. It has elements g_{ii} =  ∂ f ∂ u_{i}2 and g_{ij} = ∂ f ∂ u_{i} ⋅ ∂ f ∂ u_{j} for i,j = 1, 2, 3. The matrix product dr^{†} dr in Eq. (4.7.5) is equivalent to the dot product dr⋅dr = dr^{2} which is known as the first fundamental form or line element
dr^{2} = g_{11} du_{1}^{2} + g_{22} du_{2}^{2} + g_{33} du_{3}^{2} + 2 g_{12} du_{1} du_{2} + 2 g_{13} du_{1} du_{3} + 2 g_{23} du_{2} du_{3}  (4.7.6) 
Roughly speaking, the metric tensor g is a function which tells
how to compute the distance between any two points in a given space.
Its components can be viewed as multiplication factors which must be
placed in front of the differential displacements du_{i}
in the generalized Pythagorean theorem of Eq. (4.7.6). If u
has orthogonal Cartesian coordinates as components, i.e. r =
u, then g_{ij} =
δ_{ij} and Eq. (4.7.6) reproduces the usual
form of the Pythagorean theorem[L15Stover, C.; Weisstein, E. W. “Metric
Tensor.”
From MathWorld  A Wolfram Web Resource.]
dr^{2} = dx^{2} + dy^{2} + dz^{2}  (4.7.7) 
The volume element dV = dx dy dz can be calculated as the absolute value of the mixed product of the three tangent vectors ∂ f ∂ u_{j} (j = 1, 2, 3):
dV =  ∂ f ∂ u_{1} ⋅ ( ∂ f ∂ u_{2} × ∂ f ∂ u_{3} ) du_{1} du_{2} du_{3} = det (J) du_{1} du_{2} du_{3} = √det (g) du_{1} du_{2} du_{3}  (4.7.8) 
Now we can find the probability density function of u,
ρ(u).
Since r = f(u) establishes a
onetoone correspondence between the coordinate systems (x,
y, z) and (u_{1}, u_{2},
u_{3}), it has an inverse function such that
u = f^{−1}(r) .
In other words, the inverse function f^{−1}
(not to be confused with the reciprocal of f) is a function
that “reverses” the function f, i.e. if the function
f applied to an input u gives a result of r,
then applying its inverse function f^{−1} to
r gives the result u, and viceversa.[L16Wikipedia contributors, "Inverse function"
Wikipedia, The Free Encyclopedia.]
If ρ(r) dxdydz is the charge contained in the volume element dxdydz , applying the usubstitution yields
ρ(r) dxdydz = ρ(f(u)) det (J_{f}(u)) du_{1}du_{2}du_{3} = ρ(u) du_{1}du_{2}du_{3}  (4.7.9) 
where
ρ(u) = ρ(f(u))
det (J_{f}(u))
 (4.7.10) 
is the probability density function of the random vector u, and viceversa
ρ(r) = ρ(f^{−1}(r))
det (J_{f−1}(r))
 (4.7.11) 
is the probability density function of the random vector r.
Finally, we consider a few examples of transformations from the Cartesian coordinate system to curvilinear coordinates, i.e. to coordinate systems for Euclidean space in which the coordinate lines may be curved. Commonly used curvilinear coordinate systems include: rectangular, spherical, and cylindrical coordinate systems. These coordinates may be derived from a set of Cartesian coordinates by using a transformation that is locally invertible (a onetoone map) at each point. This means that one can convert a point given in a Cartesian coordinate system to its curvilinear coordinates and back.
r = f(u)  u = f^{−1}(r)  

coordinates 
r =

u =


Jacobian matrix 
J_{f}(u) =
∂ (x,
y, z)
∂
(R, φ, z)
=

J_{f−1}(r) =
∂ (R, φ, z)
∂ (x, y, z)
=


Jacobian (determinant) 
 ∂ (x, y, z) ∂ (R, φ, z)  = R   ∂ (R, φ, z) ∂ (x, y, z)  = 1 R  
metric tensor 
J_{f}(u)^{†}
J_{f}(u) =

J_{f−1}(r)^{†}
J_{f−1}(r) =


line element 
dr = J_{f}(u) du =
dr^{2} = dR^{2} + R^{2} dφ^{2} + dz^{2} 
du = J_{f−1}(r)
dr =
1
R^{2}
du^{2} = dx^{2} + 1 R^{2} dy^{2} + dz^{2} 

volume element 
dx dy dz = det (J_{f}(u)) dR dφ dz = R dR dφ dz  dR dφ dz = det (J_{f−1}(r)) dx dy dz = 1 R dx dy dz 
Electrostatics deals with the circumstance in which nothing depends on the time, i.e. the static case. All charges (and their density) are permanently fixed in space. The basic equation of electrostatics is Gauss' law, which in its differential form
∇ ⋅ E = 4 π ρ(r)  (E.1) 
states that the divergence of the electric field E is
proportional to the total electric charge density ρ(r).
Eq. (E.1) is written in the natural system of Hartree atomic units. In the SI
the right hand side of eq. (E.1) is multiplied by the Coulomb's constant,
k_{e}.[1(a) Mohr, P. J.;
Newell, D. B.; Taylor, B. N.
Reviews of Modern Physics 2016, 88, 035009, 73 pages.
(b) Mohr, P. J.; Newell, D. B.;
Taylor, B. N.
J. Phys. Chem. Ref. Data 2016, 45,
043102, 74 pages., N2Definition
of Hartree atomic units.]
From this fundamental equation connecting the charge density with the electric
field, the potential Φ of the field can be derived. Introducing
the potential Φ in such a way its gradient equals
−E,
E = − ∇ Φ  (E.2) 
eq. (E.1) can be rewritten in the form of a Poisson equation[N4aPoisson's equation.]
∇^{2} Φ = − 4 π ρ(r)  (E.3) 
The nabla symbol ∇ in eqs (E.13) indicates a differential vector operator whose components are partial derivative operators. It is most commonly used as a convenient mathematical notation to simplify expressions for the gradient (as in eq. E.2), divergence (as in eq. E.1), curl, directional derivative, Laplacian (as in eq. E.3), Jacobian and tensor derivative.[N4The vector differential operator ∇.]
3.1.  The electric potential. A solution of Poisson's equation[N4aSolution of Poisson's equation.] is given by eq. (E.4)
Φ(r) = ∫ ℝ^{3} ρ(r') r − r' d^{3}r'  (E.4) 
which defines the electric potential for a continous charge distribution of density ρ(r). In the case of a discrete charge distribution, where the charges are characterized by a strength q_{i} and are located at certain positions r_{i}, the charge density is given by eq. (1.2). Therefore, inserting eq. (1.2) into eq. (E.4) yields
Φ(r) = N ∑ i=1 q_{i} ∫ ℝ^{3} δ(r' − r_{i}) r − r' d^{3}r' = N ∑ i=1 q_{i} r − r_{i} = N ∑ i=1 Φ_{i}(r)  (E.5) 
where the second equality holds because of the sifting property of the Dirac delta function.[N3The Dirac Delta Function.]
The SI unit of electric potential is the volt (V, in honor of Alessandro Volta), which corresponds to J ⋅ C^{−1} and in terms of SI base units is kg ⋅ m^{2} ⋅ s^{−3} ⋅ A^{−1}. The atomic unit of electric potential (E_{h}/e) has a value of 27.111 386 02(17) V.
3.2.  The electric field. The electric field is a vector quantity generated by a distribution of electric charges, as described by Gauss' law, either in its differential form (E.1) or in its integral form.[N7Gauss' flux theorem.] Using the definition of eq. (E.2), an expression of the electric field, E(r), is derived for both a continous charge distribution
E(r) = ∫ ℝ^{3} ρ(r') r − r' r − r'^{3} d^{3}r'  (E.6) 
and a discrete charge distribution
E(r) = N ∑ i=1 q_{i} r − r_{i}^{ } r − r_{i}^{3} = N ∑ i=1 E_{i}(r)  (E.7) 
The SI units of the electric field are newtons per coulomb (N ⋅ C^{−1}) or, equivalently, volts per metre (V ⋅ m^{−1}), which in terms of SI base units are kg ⋅ m ⋅ s^{−3} ⋅ A^{−1}. The atomic unit of electric field (E_{h}/e a_{0}) has a value of 5.142 206 707(32) × 10^{11} V ⋅ m^{−1}.
3.3.  The electric field gradient. The electric field gradient is a tensor quantity generated by a distribution of electric charges. Indeed, a given charge distribution generates an electric potential Φ according to equations (E.4) and (E.5). The gradient of this potential is the negative of the electric field generated, as defined by eq. (E.2). The first derivatives of the field, or the second derivatives of the potential, yields the electric field gradient (EFG) or the hessian of the electric potential.
H(r) = − ∇
E^{†}(r) = ∇
∇^{†} Φ(r) =
 (E.8) 
H_{αβ}(r) = − ∂ E_{β}(r) ∂ r_{α} = ∂^{2} Φ(r) ∂ r_{α} ∂ r_{β} = N ∑ i=1 q_{i} [ 3 (r_{α} − r_{i,α}) (r_{β} − r_{i,β}) r − r_{i}^{5} − δ_{αβ} r − r_{i}^{3} ]  (E.9) 
Tr H(r) = ∇^{2} Φ(r) = 0  (E.10) 
The SI units of the electric field gradient are V ⋅ m^{−2}, which in terms of SI base units are kg ⋅ s^{−3} ⋅ A^{−1}. The atomic unit of electric field gradient (E_{h}/e a_{0}^{2}) has a value of 9.717 362 356(60) × 10^{21} V ⋅ m^{−2}.
3.4.  Coulomb's law. The force exerted by one point charge on another acts along the line between the charges. It varies inversely as the square of the distance separating the charges and is proportional to the product of the charges. The force is repulsive if the charges have the same sign and attractive if the charges have opposite signs. If q_{1} is at position r_{1} and q_{2} is at r_{2}, the force F_{12} exerted by q_{1} on q_{2} is the vector
F_{12} = k_{e} q_{1} q_{2} r_{12}^{3} r_{12}  (E.11) 
where r_{12} = r_{2} −
r_{1} is the vector pointing from q_{1}
to q_{2}, and k_{e} is an experimentally
determined positive constant called the Coulomb constant.[1(a) Mohr, P. J.;
Newell, D. B.; Taylor, B. N.
Reviews of Modern Physics 2016, 88, 035009, 73 pages.
(b) Mohr, P. J.; Newell, D. B.;
Taylor, B. N.
J. Phys. Chem. Ref. Data 2016, 45,
043102, 74 pages., N2k_{e} in SI units
= 8.987 551 787 3681 × 10^{9} kg m^{3}
s^{−2} C^{−2}.] The
magnitude of vector (E.11), i.e. the magnitude of the electric
force exerted by a point charge q_{1} on another point charge
q_{2} a distance r_{12} away, is given by
F_{12} = k_{e} q_{1} q_{2} r_{12}^{2}  (E.12) 
In the system of Hartree atomic units the Coulomb constant is unity by definition and therefore it will be omitted hereafter. In accord with Newton’s third law,[N7Newton’s third law.] the electrostatic force F_{21} exerted by q_{2} on q_{1} is the negative of F_{12}, i.e. F_{21} = − F_{12}.
It is worth noting that, by combining eqs (E.7) and (E.11), the electrostatic force F_{12} exerted by q_{1} on q_{2} can be written as
F_{12} = q_{2} E_{1}(r_{2})  (E.13) 
where E_{1}(r_{2}) is the electric field generated by the point charge q_{1} at position r_{2}. According to eq. (E.13), the field E_{1}(r_{2}) can be viewed as the force per unit charge on q_{2} (due to all other charges). In a system of charges, each charge exerts a force given by eq. (E.13) on every other charge. The net force on any charge q is the vector sum of the individual forces exerted on that charge by all the other charges in the system.
F(r) = q E(r)  (E.14) 
where E(r), according to eq. (E.7), is the electric field generated by the system of charges at position r. This follows from the principle of superposition of forces, which states that the force acting on a point charge due to a system of point charges is simply the vector addition of the individual forces acting alone on that point charge due to each one of the charges. The resulting force vector is parallel to the electric field vector at that point, with that point charge removed. Therefore, the force F(r) on a small charge q at position r, due to a system of N discrete charges in vacuum is:
F(r) = q N ∑ i=1 r − r_{i} r − r_{i}^{3} = q N ∑ i=1 R_{i} R_{i}^{3}  (E.15) 
where q_{i} and r_{i} are the magnitude and position respectively of the ith charge, and the vector R_{i} = r − r_{i} points from charges q_{i} to q. For a continous charge distribution where ρ(r′) gives the charge per unit volume at position r′, the principle of superposition is also used. In this case an integral over the region containing the charge is equivalent to an infinite summation over infinitesimal elements of space containing a point charge dq' = ρ(r′) d^{3}r' (where d V' = d^{3}r' is an infinitesimal element of volume):
F(r) = q ∫ ℝ^{3} ρ(r') r − r' r − r'^{3} d^{3}r'  (E.16) 
There are three conditions to be fulfilled for the validity of Coulomb’s
law:[L?5Coulomb’s law
From Wikipedia, the free encyclopedia]
The SI unit of force is the newton (N), which in terms of SI base units is kg ⋅ m ⋅ s^{−2}. In atomic units the force is expressed in hartrees per Bohr radius (E_{h}/a_{0}). The atomic unit of force has a value of 8.238 723 36(10) × 10^{−8} N. The CGS unit of force is the dyne, which in in terms of CGS base units is g ⋅ cm ⋅ s^{−2}, so the SI unit of force, the newton, is equal to 10^{5} dynes.
3.5.  Electrostatic potential energy and electrostatic potential energy density. Electrostatic potential energy, U_{E}, is the energy possessed by an electrically charged object because of its position relative to other electrically charged objects.
The electrostatic potential energy, U_{E}, of one point charge q at position r in the presence of an electric field E is defined as the negative of the work W done by the electrostatic force (E.14) to bring q from infinity to position r,
U_{E} = − W_{∞→r} = − q r ∫ ∞ E(r') ⋅ d r' = q r ∫ ∞ ∇ Φ(r') ⋅ d r' = q [Φ(r) − Φ(∞)] = q Φ(r)  (E.17) 
where E is electrostatic field defined by eq. (E.2) and dr' is the displacement vector in a curve from infinity to the final position r. The third equality holds beacause Coulomb forces are conservative, so that the work they do in moving a particle between two points is independent of the taken path. In the limit as r → ∞, the electrostatic potential Φ(r) falls off like 1/r, and is consequently zero. This explains the last equality in eq. (E.17).
U_{E} = ∑ all pairs q_{i} q_{j} r_{ij} = 1 2 N ∑ i=1 q_{i} N ∑ j≠i q_{j} r_{ij} = 1 2 N ∑ i=1 q_{i} Φ_{i}  (E.18) 
U_{E}  = 1 2 (U_{E} + U_{E}) 
= 1 2 ( q_{2} q_{1} r_{12} + q_{3} q_{1} r_{13} + q_{3} q_{2} r_{23} + q_{2} q_{1} r_{12} + q_{3} q_{1} r_{13} + q_{3} q_{2} r_{23} )  
= 1 2 [ q_{1} ( q_{2} r_{12} + q_{3} r_{13} ) + q_{2} ( q_{3} r_{23} + q_{1} r_{12} ) + q_{3} ( q_{1} r_{13} + q_{2} r_{23} ) ]  
= 1 2 (q_{1} Φ_{1} + q_{2} Φ_{2} + q_{3} Φ_{3}) 
Equation (E.18) also describes the electrostatic potential energy of a continuous charge distribution. By considering that each volume element dV contains the element of charge dq = ρ dV, eq. (E.18) can be written
U_{E} = 1 2 ∫ all space ρ(r_{1}) ρ(r_{2}) r_{12} dV_{1} dV_{2} = 1 2 ∫ V_{1} ρ(r_{1}) Φ(r_{1}) dV_{1} ≡ 1 2 ∫ V ρ(r) Φ(r) dV  (E.19) 
where the factor ^{1}⁄_{2} is introduced because in the double integral over dV_{1} and dV_{2} all pairs of charge elements have been counted twice. In addition, the second equality in eq. (E.19) follows from the fact that the integral over dV_{2} is just the potential at r_{1} (cf eq. E.4):
Φ_{1} ≡ Φ(r_{1}) = ∫ V_{2} ρ(r_{2}) r_{12} dV_{2}  (E.20) 
The last identity in eq. (E.19) comes from the fact that the point r_{2} no longer appears in the second equality, so that the suffix 1 is reduntant and can be safely neglected.
Finally, it's worth to note that the electrostatic energy can be described not only in terms of the charges, but also in terms of the electric fields they produce:
U_{E} = ∫ V u_{E} dV = 1 8 π ∫ V E^{2}(r) d^{3}r  (E.21) 
This formula shows that when an electric field E(r) is present, there is located in space an energy U_{E} whose density (energy per unit volume) is proportional to the square of the electric field strength, i.e. u_{E} = E^{2}(r) ⁄ 8π.
U_{E} = − 1 8 π ∫ V Φ ∇^{2}Φ dV  (E.211) 
Φ ∇^{2}Φ  = Φ ( ∂^{2} Φ ∂ x^{2} + ∂^{2} Φ ∂ y^{2} + ∂^{2} Φ ∂ z^{2} )  
= ∂ ∂ x ( Φ ∂ Φ ∂ x ) − ( ∂ Φ ∂ x )^{2} + ∂ ∂ y ( Φ ∂ Φ ∂ y ) − ( ∂ Φ ∂ y )^{2} + ∂ ∂ z ( Φ ∂ Φ ∂ z ) − ( Φ ∂ Φ ∂ z )^{2}  (E.212)  
= ∇ ⋅ (Φ ∇Φ) − (∇Φ) ⋅ (∇Φ) 
U_{E} = 1 8 π ∫ V (∇Φ) ⋅ (∇Φ) dV − 1 8 π ∫ V ∇ ⋅ (Φ ∇Φ) dV  (E.213) 
∫ V ∇ ⋅ (Φ ∇Φ) dV = ∫ S (Φ ∇Φ) ⋅ n dS  (E.214) 
U_{E} = 1 8 π ∫ V (∇Φ) ⋅ (∇Φ) dV = 1 8 π ∫ V E ⋅ E dV  (E.215) 
The SI unit of energy is the joule (J), which in terms of SI base units is kg ⋅ m^{2} ⋅ s^{−2}. In atomic units the energy is expressed in hartree (E_{h}). The atomic unit of energy has a value of 4.359 744 17(75) × 10^{−18} J or, in more common non SI units, 27.211 385 eV and 627.509 kcal ⋅ mol^{−1}.
Definition of the Dirac delta function in three
dimensions. The definition of the Dirac delta function in one
dimension, δ(x − x_{0}), can
be found in reference [4Messiah, A.
“Quantum Mechanics”, Volume I
NorthHolland Publishing
Company, Amsterdam: 1970, pp. 179184 and 468471.].
Here its definition is extended to the threedimensional case, where the
delta function can be indicated by
δ_{3}(r −
r_{0}) or simply δ(r
− r_{0}) in terms of the two position
vectors r and r_{0}.
The easiest way to define a threedimensional delta
function is just to take the product of three onedimensional functions,
i.e. in cartesian coordinates:
δ(r − r_{0}) = δ(x − x_{0}) δ(y − y_{0}) δ(z − z_{0})  (N3.1) 
Expressions for cylindrical and spherical coordinates can be be found at
link [L3Weisstein, Eric W.
“Delta Function.”
From MathWorld − A Wolfram Web
Resource.
http://mathworld.wolfram.com/DeltaFunction.html].
The Dirac delta function δ(r −
r_{0}) is defined by the property
∫ f(r) δ(r − r_{0}) d^{3}r = f(r_{0})  (N3.2) 
for any function f(r) that is continous at the point r = r_{0}. The integral (N3.2) is over any volume that contains the point r = r_{0}. At this point δ(r − r_{0}) is singular in such a manner that
∫
δ(r − r_{0})
d^{3}r = ∞ ∫ −∞ δ(x − x_{0}) dx ∞ ∫ −∞ δ(y − y_{0}) dy ∞ ∫ −∞ δ(z − z_{0}) dz = 1  (N3.3) 
In addition, δ(r − r_{0}) is zero everywhere except at r = r_{0}, i.e.
δ(r − r_{0}) = 0 ∀ r ≠ r_{0}  (N3.4) 
The discrete analog of the Dirac delta function is the Kronecker delta function, which is usually defined on a discrete domain and takes values 0 and 1.
Eq. (N3.2) is also known as the “sifting” property of the Dirac delta function since that, when the delta function appears in a product within an integrand, it has the property of “sifting” out the value of the integrand at the point where the delta function is not null.
∇ =
∂
∂ r
=
 (N4.1) 
Mainly, it is a convenient mathematical notation that allows to write, with compact and intuitive notation, different differential operators, like the gradient, the divergence and the curl, depending on the way it is applied.
When applied to a continuous function f: ℝ^{3} → ℝ, i.e. a function f(r) which operates on a vector variable r = (x, y, z)^{†} ∈ ℝ^{3} and returns a scalar, it defines the gradient of the scalar field f = f(r):
grad f =
∇ f =
 (N4.2) 
The definition (N4.2) can be extended to a coninuous vector function f: ℝ^{3} → ℝ^{m} which takes as input the vector r ∈ ℝ^{3} and produces as output f(r) ∈ ℝ^{m}. Then the gradient of the vector field f = f(r) is the outer product or tensor product of the differential vector operator ∇^{†} with the vector f, which yields a tensor quantity, i.e. the Jacobian m×3 matrix J_{f} of f, defined and arranged as follows:
grad f = ∇^{†}⊗f =
 (N4.3) 
Similarly, if F: ℝ^{3} → ℝ^{m}×ℝ^{n} is a continuous second order tensor function which takes as input the vector r ∈ ℝ^{3} and produces as output the m×n matrix F(r) ∈ ℝ^{m}×ℝ^{n}, then the gradient of F is the Jacobian m×3n matrix J_{F} of F.
Eqs. (N4.2) and (N4.3) show that, instead of “grad”, it is worth using the operator ∇^{†}⊗ where the symbol ⊗ stands for the Kronecker product of two matrices[N6Kronecker product (definition)], and is a generalization of the outer product or tensor product of vectors.
The divergence of a continuously differentiable vector field f = (f_{x}, f_{y}, f_{z})^{†} is defined as the scalar field:
div f = ∂ f_{x} ∂ x + ∂ f_{y} ∂ y + ∂ f_{z} ∂ z = ∇ ⋅ f  (N4.4) 
where the last equality, ∇ ⋅ f, is the
common notation for the divergence, and is a convenient mnemonic,
reminescent of the dot product (also called scalar product
or, less frequently, inner product) of two vectors:
take the components of the ∇ operator, apply
them to the corresponding components of f, and sum the results.[L11Wikipedia contributors, "Divergence"
Wikipedia, The Free Encyclopedia.]
Because applying an operator is different from multiplying the components,
this is considered an abuse of notation.[L12Wikipedia contributors, "Abuse of notation"
Wikipedia, The Free Encyclopedia., L13Tai, Chen "A survey of the improper use of ∇
in vector analysis" (1994).]
The dot product ∇ ⋅ f is equivalent to a
matrix multiplication ∇^{†} f,
provided that ∇ is represented as a 3×1 column
matrix (which makes ∇^{†} a 1×3
row matrix) and f as a 3×1 column matrix:
div f = ∇^{†} f =
 (N4.5) 
It is worth noting (compare, for example, Eqs (N4.4) and (N4.3)) that the inner product for Euclidean vector spaces is the trace of the outer product.
Notation (N4.5) can also be used to define the divergence of a
continuously differentiable secondorder tensor field F, which
turns out to be a firstorder tensor field:[L11Wikipedia contributors, "Divergence"
Wikipedia, The Free Encyclopedia.]
div F = ∇^{†} F =
 (N4.6) 
where the elements of the firstorder tensor field (a vector field) in the last equality are the divergence of the 3×1 column matrices F_{i} = (F_{xi}, F_{yi}, F_{zi})^{†} (i = x, y, z) extracted from matrix F.
The rotation (curl f or rot f) of the vector field f is a vector, which is symbolically expressed by the cross product of ∇ and the vector f:
curl f = ∇ × f
= det
 (N4.7) 
The Laplacian (div grad f) of the scalar field f is a scalar, which is symbolically expressed by the scalar multiplication of ∇^{2} and the scalar field f:
div grad f = ∇^{†} ∇ f = ∇ ⋅ ∇ f = ∇^{2} f = ∂^{2} f ∂ x^{2} + ∂^{2} f ∂ y^{2} + ∂^{2} f ∂ z^{2}  (N4.8) 
A matrix which has only one column is called a column matrix. The transpose of a column matrix is another matrix which has only one row and is called a row matrix. Since a matrix is defined as a bidimensional array whose elements have exactly two indices, a column matrix has dimensions n×1, and a row matrix has dimensions 1×n. In spite of this, they are often indicated as vectors, which instead are defined as monodimensional arrays of elements with exactly one index. Programming languages like Fortran that support multidimensional arrays typically have a native columnmajor or canonical storage order for these arrays, that means that consecutive elements of the columns of the arrays are contigous in memory. This linear storage of multidimensional arrays or tensors complies with tensors being represented by a column matrix.
A⊗B ≔

The Poisson's equation is a partial differential equation (PDE) of the form
∇^{2} u(r = f(r)  (N6.1) 
which describes the response u(r) of a physical system to an arbitrary source f(r). The general solution of eq. (N6.1) is given by the sum of any harmonic function [which satisfies the homogeneous Laplace equation, ∇^{2} u(r) = 0] and a particular solution of the nonhomogeneous Poisson equation (N6.1). Let's first consider the nonhomogeneous problem associated with the presence of a point source located at r', mathematically represented by the Dirac delta function δ(r − r'):
∇^{2} G(r, r') = δ(r − r')  (N6.2) 
where the Green's function G(r, r') describes the response of the system at the point r to a point source located at r'.
If G is the fundamental solution of eq. (N6.2), then the convolution
G * f = ∫ ℝ^{3} G(r, r') f(r') d^{3}r'  (N6.3) 
will be a solution of eq. (N6.1). Indeed, because the Laplacian ∇^{2} is a linear differential operator, it follows that
∇^{2}(G * f ) =
(∇^{2} G) * f =
δ * f =
∫ ℝ^{3} δ(r − r') f(r') d^{3}r' = f(r)  (N6.4) 
where the first equality comes from linearity of the Laplacian, the second equality comes from eq. (N6.2), the third equality follows from the definition of convolution, and the last equality is a consequence of the fact that the Dirac delta function is an identity element for convolution, as stated by eq. (N6.2). Therefore, if G is the fundamental solution, the convolution u = G * f is one solution of ∇^{2} u(r) = f(r). This does not mean that it is the only solution. Several solutions for different initial conditions can be found.
The definition (N6.2) of the fundamental solution implies that, if the Laplacian of G is integrated over any volume that encloses the source point, then, according to eq. (N6.3),
∭ V ∇^{2} G dV = 1  (N6.5) 
The Laplace equation is unchanged under a rotation of coordinates, and hence we can expect that a fundamental solution may be obtained among solutions that only depend upon the distance r from the source point. If we choose the volume to be a ball of radius a around the source point, then Gauss' divergence theorem[N6aThe divergence theorem.] implies that
∭ V ∇^{2} G dV = ∬ S d G d r d S = 4 π a^{2} d G d r _{r=a} = 1  (N6.6) 
Rearranging the last equality of eq. (N6.6) yields
d G d r _{r=a} = 1 4 π a^{2}  (N6.7) 
and hence
G(r, r') = − 1 4 π a = − 1 4 π ⋅ 1 r − r'  (N6.8) 
as a^{2} = (x − x')^{2} + (y − y')^{2} + (z − z')^{2} = r − r'^{2} is the equation of a sphere of radius a centered at the source point r'. Finally, inserting eq. (N6.8) into eq. (N6.3) gives the solution of Poisson's equation (N6.1):
u(r) = − 1 4 π ∫ ℝ^{3} f(r') r − r' d^{3}r'  (N6.9) 
Eq. (N6.9) is a solution of Poisson's equation of electrostatics (E.3), provided that u(r) = Φ(r) and f(r) = − 4 π ρ(r).
The divergence theorem states that the surface integral of the normal component of an arbitrary vector F over a closed surface S is equal to the volume integral of the divergence of the vector over the volume V interior to the surface S:
∬ S (F ⋅ n) dS = ∭ V (∇ ⋅ F) dV  (N7.1) 
where dS may be used as a shorthand for n dS. The lefthand side of the equation represents the total flow across the closed boundary S, and the righthand side represents the total of the sources in the volume V. In other words, the divergence theorem is a mathematical statement of the physical fact that, in the absence of the creation or destruction of matter, the density within a region of space can change only by having it flow into or away from the region through its boundary. In the special case where the vector F is defined as the gradient of a scalar G, F = ∇ G, and the volume V is that is that of a ball of radius a, the divergence theorem yields eq. (N6.5).
Newton's third law: forces always occur in equal and opposite pairs. If object A exerts a force F_{A,B} on object B, an equal but opposite force F_{B,A} is exerted by object B on object A. Thus,
F_{B,A} = − F_{A,B}  (N8.1) 
This is also known as the Action  Reaction Principle: "For every action there is always an opposite and equal reaction".