# Cartesian Multipole Moment Tensors and Polytensors.

## 1. - Unabridged Moments

The unabridged cartesian moment tensor of order n of a charge distribution ρ(r) is defined by:

 m(n) = ∫ρ(r) r(n) d3r ≡ ⟨r(n)⟩ (1.1)

where the integral is over the volume containing the charge distribution and the moment tensor operator r(n) denotes the tensor product of vector r by itself n times:

 r(n) = r⊗r⊗…   (n times) (1.2)

For this reason, r(n) is also called the n-th tensor power of the vector r and is some times indicated as rn. The symbol ⊗ denotes the Kronecker product of two matrices of arbitrary size.[1Kronecker product (definition), 2Zhang, H.; Ding, F.
Journal of Applied Mathematics 2013, Article ID 296185, 8 pages.
]

The elements of the tensor r(n) are all the products of n factors rα1rα2rα3…rαn (usually indicated by the index sets α1α2α3…αn), where each subscript αi denotes Cartesian axes x, y, and z. The number of the elements of the tensor r(n) is the number of n-tuples (i.e., ordered lists of length n) of a set of three elements x, y, and z and is equal to 3n. The element ordering of r(n) is usually chosen to be canonical,[3Applequist, J.
J. Math. Phys. 1983, 24, 736-741.
] which means that the first index is changing fastest (i.e., on proceeding through the array, α1 varies through the values 1,2,3 more rapidly than α2, which varies more rapidly than α3, etc.). The canonical order is the order in which array elements are usually stored in a computer memory. On the other hand, the elements are said to be in anticanonical order when the first index is changing slowest. The elements of the tensor r(n) comprise a column matrix.[4Definition of column matrix] To emphasize this fact, the tensor r(n) of order n may be represented as a tensor of subdivided order (n,0) or (0,n), which corresponds to a column or row matrix, respectively, i.e. r(n)r(n,0)r(0,n). In general, any tensor T(m+n) of order m + n may be represented as a tensor of subdivided order (m,n), also called a tensor of type (m,n), T(m,n), in which m and n are the order indices. The array of components, T(m,n)α1α2α3…αmβ1β2β3…βn, comprise a rectangular matrix whose rows are indexed by the set α1α2α3…αm and whose columns are indexed by the set β1β2β3…βn.[3Applequist, J.
J. Math. Phys. 1983, 24, 736-741.
]

For n = 0, the zero-th order moment operator is the unit operator, r(0) = 1 (a scalar), and for n = 1, the first (order) moment operator is the position vector:

r(1)r =
 x y z
(1.3)

For n = 2, the second (order) moment operator r(2) is the dyadic (tensor) product rr:[5Dyadic product of two vectors of the same dimension]

r(2) = rr =
 x r y r z r
=
 xx xy xz yx yy yz zx zy zz
= vec
 xx yx zx xy yy zy xz yz zz
≡ vec(Q) =
(1.4)
= vec Q = vec
 xx xy xz yx yy yz zx zy zz
= vec
 x r† y r† z r†
= vec(rr) =
 xx yx zx xy yy zy xz yz zz

The first line of Eq. (1.4) shows that the Kronecker product of the column matrix r(1,0)r by itself yields, by definition,[1Kronecker product] a column matrix r(2,0) whose elements are in anticanonical order. Such a column matrix, which represents a tensor of subdivided order (2,0), can be viewed as the vectorization of the 3×3 quadrupole matrix operator Q, which represents a tensor of type (1,1). The vectorization is accomplished by using the “vec” operator, which stacks the columns of a matrix one underneath the other to form a single vector.[6Magnus, J. R.; Neudecker, H.
The Annals of Statistics 1979, 7, 381-394.
, 7Henderson, H. V.; Searle, S. R.
Linear and Multilinear Algebra 1981, 9, 271-288.
] Since matrix Q is symmetric, i.e. Q = Q, the equalities of the second line of Eq. (1.4) follow straightforwardly. In other words, Eq. (1.4) states that the 3×3 matrix rr, obtained by multiplying each element of r by each element of r or, equivalently, by the Kronecker product rr, gives directly, after vectorization, a tensor of type (2,0) with elements in canonical order. The index notation (rr)ij = rirj = rjri = (rr)ji, where both i and j run over 1, 2 and 3 or x, y and z, emphasizes the fact that r(2) is a symmetric tensor. However, because the tensor is symmetric and rr = (rr), the component indexes can be safely exchanged to move from canonical to anticanonical order, and viceversa.

Similarly, for n = 3, the third (order) moment operator is a 3×3×3 tridimensional array, which can be written as:

r(3) = r(2)r = vec(rr)⊗r =
 xx r† yx r† zx r† xy r† yy r† zy r† xz r† yz r† zz r†
=
 xxx xxy xxz yxx yxy yxz zxx zxy zxz xyx xyy xyz yyx yyy yyz zyx zyy zyz xzx xzy xzz yzx yzy yzz zzx zzy zzz
= r(2,1)
(1.5)
 xxx xyx xzx xxy xyy xzy xxz xyz xzz yxx yyx yzx yxy yyy yzy yxz yyz yzz zxx zyx zzx zxy zyy zzy zxz zyz zzz
= r(1,2)

It is represented by either a 9×3 matrix r(2,1) or, equivalently, by a 3×9 matrix r(1,2), that are tensors of type (2,1) and (1,2), respectively. Both matrices are stored columnwise in canonical order, so that they are implicitly vectorized as a r(3,0) column matrix. In index notation, the tensor is indicated by   [r(3)]ijk = rirjrk,   where i, j and k run over 1, 2, and 3, i.e. over x, y and z. Due to the symmetry of the tensor, the ordering of the components can be changed from canonical to anticanonical mode by simply reversing the indices.

In general, the procedure outlined in Eqs. (1.4) and (1.5) can be used to build up the n-th moment operator:

 r(n) = r(n−1)⊗r† (1.6)

where r(n) and r(n−1) are tensors of type (n,0) and (n−1,0), respectively.

A computer code for the generation of the component labels in canonical order of the n-th order moment operator r(n) may require n nested do-loops. The Program 1 provides an example of Fortran 77 code for n = 3.

 ```Character*1 R(3) /'x', 'y', 'z'/ Character*3 S(3**3) ijk = 0 Do k=1,3 Do j=1,3 Do i=1,3 ijk = ijk + 1 S(ijk)(1:1) = R(i) S(ijk)(2:2) = R(j) S(ijk)(3:3) = R(k) end Do end Do end Do Write(*,'(3(7(a3,", "),/),5(a3,", "),a3)') (S(m),m=1,27) ```

Generalization to any order n is provided by the Fortran 90 code of Program 2, which implements the Nested Summation Symbol (NSS) operator:

 ∑ n (j = i, f) := f1 ∑ j1=i1 f2 ∑ j2=i2 … fn ∑ jn=in (1.7)

introduced by Carbó and Besalú.[9Carbó, R.; Besalú, E.
J. Math. Chem. 1993, 13, 331-342.
, 10Carbó, R.; Besalú, E.
J. Math. Chem. 1995, 18, 37-72.
] The implied n-dimensional vectors of Eq. (1.8) are j = (j1, j2, …, jn), i = (i1, i2, …, in), and f = (f1, f2, …, fn). The index n is called the dimension of the NSS. Program 2 is a translation of the NSS into a generalized nested do-loop (GNDL) structure, which is independent of the dimension of the involved nested sums.

 ```! ------------------------------------------------------------------ Module NSS ! ------------------------------------------------------------------ implicit none private public CTLbl, PrtOut contains ! ------------------------------------------------------------------ Subroutine CTLbl (j,S) ! ------------------------------------------------------------------ Integer, intent(out) :: j(:) Character(len=1), intent(out) :: S(:,:) Character(len=1), dimension(3), parameter :: R = (/'x', 'y', 'z'/) Integer :: n, k, num n = size(j) ! initial NSS parameter values do k=1,n j(k) = 1 end do num = 0 ! count the components (num = 3^n) ! Start GNDL procedure k = n ! the innermost loop corresponds to the index k = n do while (k > 0) if (j(k) > 3) then ! index out of range j(k) = 1 ! reset this level k = k - 1 ! previous level will be activated, if possible else ! calculational kernel num = num + 1 do k=1,n S(n-k+1,num) = R(j(k)) ! canonical order end do ! end of calculational kernel k = n end if if (k > 0) then j(k) = j(k) + 1 ! step increment end if end do ! End of GNDL procedure End Subroutine CTLbl ! ------------------------------------------------------------------ Subroutine PrtOut (n,k,S) ! ------------------------------------------------------------------ Integer, intent(in) :: n, k Character(len=1), intent(in) :: S(n,k) Integer, parameter :: rowlength = 72 Integer :: l, m, fulllines, lastlinelength Character(len=80) :: fmt m = (rowlength+1)/(n+2) fulllines = k/m if (fulllines*m == k) then fulllines = fulllines - 1 lastlinelength = m - 1 if (fulllines >= 1) then write(fmt,'(1h(,i6,1h(,i2,1h(,i2,12ha1,", "),/),,i2,1h(,i2,& 9ha1,", "),,i2,3ha1))') fulllines,m,n,lastlinelength,n,n else write(fmt,'(1h(,i2,1h(,i2,9ha1,", "),,i2,3ha1))') & lastlinelength,n,n end if else lastlinelength = k - fulllines*m - 1 if (fulllines >= 1) then if (lastlinelength >= 1) then write(fmt,'(1h(,i6,1h(,i2,1h(,i2,12ha1,", "),/),,i2, & 1h(,i2,9ha1,", "),,i2,3ha1))') fulllines,m,n, & lastlinelength,n,n else write(fmt,'(1h(,i6,1h(,i2,1h(,i2,12ha1,", "),/),,i2, & 3ha1))') fulllines,m,n,n end if else write(fmt,'(1h(,i2,1h(,i2,9ha1,", "),,i2,3ha1))') & lastlinelength,n,n end if end if write(*,fmt) ((S(l,m),l=1,n),m=1,k) End Subroutine PrtOut End Module NSS ```

Cartesian moment tensors of order n are also called 2n-pole moments, i.e. mono-, di-, quadru-, octa-, and hexadeca-pole for n = 0, 1, 2, 3, 4, respectively.

Applequist has set forth a useful organization for electrostatic properties,[3Applequist, J.
J. Math. Phys. 1983, 24, 736-741.
] which is followed in PAMoC. Applequist defines a Cartesian polytensor of the first order as a sequence of Cartesian tensors, each arranged in a column. The first-degree polytensor is a stacked list of the (column) moment tensors in increasing order, and will be designated M

M =
 m(0) m(1) m(2) m(3) m(4) ⋮
.
(1.8)

The sequence of component tensors, m(i) is of indefinite length, but PAMoC for practical purposes truncates the sequence at the fourth order. Since the number of components of the cartesian moment tensor m(n) of order n is 3n, the number of components of a polytensor which contains a sequence of cartesian moment tensors up to the n-th order is n i=0 3i = (3n+1 − 1) / 2. It is useful to relate the indeces α1α2αn of an element entry of a tensor of order n to a single-index reference kn which represents the location in memory after the element of the tensor. For a complete tensor in canonical order it is

 kn = α1 + n ∑ m=2 3m−1 (αm − 1)      ( ∀ α1α2…αn = 1, 2, 3;    kn ∈ [1, 3n] ) (1.9)

and for the corresponding polytensor it is    K = kn + (3n − 1)/2.

## 2. - Symmetric and Compressed Tensors

Each component of the moment tensor operator r(n) is invariant with respect to permutations of its suffixes, i.e. the tensor is totally symmetric. Then, suffixes can be ordered according to their type, i.e. mx1xnxy1ynyz1znz, where nx, ny, and nz, are the number of times x, y, and z occur in the component index α1α2…αn. The ni are called degree indices and satisfy nx + ny + nz = n. The following concise notation can be used

 mnxnynz ≡ ⟨xnxynyznz⟩      ∀ nx, ny, nz = 0,1,2,…,n   with   nx + ny + nz = n (2.1)

where only the three degree indices are implied.

The number of distinct components of a symmetric tensor of order n is ½(n+1)(n+2). A symmetric tensor which comprises only its distinct components is said to be compressed and is usually arranged as a column vector with components in canonical order. Clearly, compression has significance only if n ≥ 2. The canonical array of index sets of the compressed tensor is derived from that of the complete tensor by deleting any index set which does not satisfy the condition   α1≥α2≥…≥αn−1≥αn. The corresponding Fortran 90 code is reported in Program 3.

 ```! ------------------------------------------------------------------ program GNDL ! ------------------------------------------------------------------ ! ! Nested Summation Symbol (NSS): sum over n (j = i, f) ! Translation of the NSS into a Generalized Nested Do-Loop (GNDL) ! ! ------------------------------------------------------------------ Use NSS Integer, allocatable, dimension(:) :: j Character(len=1), allocatable, dimension(:,:) :: S Character(len=3) :: t Integer :: k, l, m, n Logical :: ok if (IArgC() /= 1) then write(*,'(/," Enter the order n of the moment operator!",/)') stop else call GetArg (1,t) read(t,*) n if (n <= 0) then ! the order is out of range write(*,'(" Error: order",i3," out of range.")') n stop end if end If allocate ( j(n), S(n,3**n) ) call CTLbl (j,S) ! Output the complete tensor write(*,'(/," Full Cartesian Tensor of order n =",i2,/, & " with components in canonical order",/)') n call PrtOut (n,3**n,S) ! Compressed tensor k = 0 do m=1,3**n ok = .true. do l=1,n-1 ok = ok .and. S(l,m) >= S(l+1,m) end do if (ok) then k = k + 1 do l=1,n S(l,k) = S(l,m) end do end if end do ! Output the compressed tensor write(*,'(/," Compressed Cartesian Tensor of order n =",i2,/, & " with components in canonical order",/)') n call PrtOut (n,(n+1)*(n+2)/2,S) deallocate ( j, S ) End Program ```

Since the number of components of the symmetric compressed tensor m(n) of order n is (n + 1)(n + 2)/2, the number of components of a polytensor which contains a sequence of symmetric compressed tensors up to the n-th order is    n i=0 (i + 1)(i + 2)/2 = (n + 1)(n + 2)(n + 3)/6 .   The degree indeces nxnynz of a compressed canonical totally symmetric Cartesian tensor of order n = nx + ny + nz can be related to a single index reference jn

 jn = 1 + nz + (ny + nz)(ny + nz + 1)/2      ( ∀ ny, nz = 1, 2, 3;    jn ∈ [1, (n + 1)(n + 2)/2] ) (2.2)

which represents the location in memory after the element of the tensor. Its position in the corresponding polytensor is    J = jn + n(n + 1) / 2.

The number of components of complete and compressed Cartesian tensors and polytensors as a function of their order n is reported in Table 2.1.

Table 2.1 − Number of components of complete and compressed Cartesian tensors and polytensors as a function of their order n.
n complete tensor
3n
complete polytensor
(3n+1 − 1)/2
compressed tensor
(n + 1)(n + 2)/2
compressed polytensor
(n + 1)(n + 2)(n + 3)/6
0 1 1 1 1
1 3 4 3 4
2 9 13 6 10
3 27 40 10 20
4 81 121 15 35

The number of times that each of the (n+1)(n+2)/2 distinct components of a totally symmetric Cartesian tensor, m(n), appears in the complete set of 3n components is given by the number of permutations of n elements with repetitions of nx, ny and nz elements, i.e.

 gj(n) = n! / nx! ny! nz!             ∀ j ∈ [1, (n+1)(n+2)/2]   and   n = nx + ny + nz (2.3)

where j spans the compressed canonical array and ni is the number of times i appears in the component index of the j-th element. Eq. (2.3) defines the elements of the diagonal matrix g(n). It is easily verified that Tr [g(n)] = 3n. The values of gj(n) are reported in Table 2.2 for the compressed canonical tensors of order 0 through 4.

The component labels of complete Cartesian tensors of order 0 through 4 are reported in Table 2.2, using canonical ordering with both component indeces, α1α2αn, and degree indeces, nxnynz. The components are numbered by the tensor single-index kn, defined by eq. (1.9) and by the associated polytensor single-index K. The single-indeces jn and J of the canonical compressed tensors are also reported, together with their multiplicities, defined by eq. (2.3). In addition, the ordering of compressed tensors adopted in PAMoC is shown. PAMoC reports the components of the quadrupole moment in the order they are obtained by vectorizing the upper triangular part of the symmetric quadrupole matrix given by Eq. (1.4). For the component ordering of the octupole and hexadecapole moments, PAMoC uses the same convention adopted by the electronic structure modeling package GAUSSIAN-09.[11Gaussian 09, Revision E.01, Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian, Inc., Wallingford CT, 2009.]

Table 2.2 − Component indeces, α1α2αn, and degree indeces, nxnynz, of the elements of complete and compressed Cartesian tensors of order 0 through 4, numbered by tensor and polytensor single-indeces.
complete tensor
(canonical order)
compressed tensor (canonical order,
α1≥α2≥…≥αn−1≥αn)
compressed tensor
(PAMoC order)
tensor order
n
polytensor
single index,
K
tensor
single index,
kn
component
index
α1α2αn
degree index
nxnynz
polytensor
single index,
J
tensor
single index,
jn
gin(n)   polytensor
single index,
I
tensor
single index,
in
gin(n)
0 11 000 11 1 11 1
1 21x100 21 1 21 1
32y010 32 1 32 1
43z001 43 1 43 1
2 51xx200 51 1 51 1
62yx110 62 2
73zx101 73 2
84xy110 62 2
95yy020 84 1 73 1
106zy011 95 2
117xz101 84 2
128yz011 95 2
139zz002 106 1 106 1
3 14 1xxx300 111 1 11 1 1
15 2yxx210 122 3
16 3zxx201 133 3
17 4xyx210
18 5yyx120 144 3
19 6zyx111 155 6
20 7xzx201
21 8yzx111
22 9zzx102 166 3
2310xxy210 15 5 3
2411yxy120
2512zxy111
2613xyy120 14 4 3
2714yyy030 177 1 12 2 1
2815zyy021 188 3
2916xzy111
3017yzy021
3118zzy012 199 3
3219xxz201 16 6 3
3320yxz111
3421zxz102
3522xyz111 2010 6
3623yyz021 19 9 3
3724zyz012
3825xzz102 17 7 3
3926yzz012 18 8 3
4027zzz003 20101 133 1
4 41 1xxxx400 21 1 1 21 1 1
42 2yxxx310 22 2 4
43 3zxxx301 24 3 4
44 4xyxx310
45 5yyxx220 24 4 6
46 6zyxx211 25 512
47 7xzxx301
48 8yzxx211
49 9zzxx202 26 6 6
5010xxyx310
5111yxyx220
5212zxyx211
5313xyyx220
5414yyyx130 27 7 4 26 6 4
5515zyyx121 28 812
5616xzyx211
5717yzyx121
5818zzyx112 29 912
5919xxzx301
6020yxzx211
6121zxzx202
6222xyzx211
6323yyzx121
6424zyzx112
6525xzzx202
6626yzzx112
6727zzzx103 3010 4 28 8 4
6828xxxy310 24 4 4
6929yxxy220
7030zxxy211
7131xyxy220
7232yyxy130
7333zyxy121
7434xzxy211
7535yzxy121
7636zzxy112 351512
7737xxyy220 3010 6
7838yxyy130
7939zxyy121
8040xyyy130
8141yyyy040 3111 1 22 2 1
8242zyyy031 3212 4
8343xzyy121
8444yzyy031
8545zzyy022 3313 6
8646yxzy121
8747zxzy112
8848xxzy211
8949yyzy031
9050zyzy022
9151xyzy121
9252yzzy022
9353zzzy013 3414 4 29 9 4
9454xzzy112
9555yxxz211
9656zxxz202
9757xxxz301 25 5 4
9858yyxz121 341412
9959zyxz112
10060xyxz211
10161xzxz202
10262yzxz112
10363zzxz103
10464xxyz211 331312
10565yxyz121
10666zxyz112
10767xyyz121
10868yyyz031 27 7 4
10969zyyz022
11070xzyz112
11171yzyz022
11272zzyz013
11373xxzz202 3111 6
11474yxzz112
11575zxzz103
11676xyzz112
11777yyzz022 3212 6
11878zyzz013
11979xzzz103
12080yzzz013
12181zzzz004 3515 1 23 3 1

From an algebraic point of view, a (n+1)(n+2)/2 × 3n matrix C can be defined to transform a complete canonical tensor, m(n), into a compressed canonical tensor, m(n):

 m(n) = C m(n) (2.4)

with the inverse transformation:[20 The product C C is not positive definite and cannot be inverted.
]

 m(n) = C† g(n) m(n) (2.5)

The compression matrix C is defined in such a way:

 C C† = [g(n)]−1     and     (C C†)−1 = g(n) (2.6)

i.e. let the row-suffix k of a generic element Ckl specify the canonical sequence number of the compressed-tensor component, whereas the column-suffix l is the canonical sequence number of a complete-tensor component. Then  Ckl = gk−1  for all values of l which identify any of the gk components of the complete tensor equal by symmetry to the k-th component of the compressed tensor. Otherwise, Ckl = 0. However, it is worth noting that PAMoC never uses matrix C, explicitly.

## 3. - Translation of Unabridged Moments

Consider the transformation of the unabridged cartesian moment tensor operator of order 1 (position vector) under translation of the axis system given by

 r = r − R (3.1)

In eq. (3.1), r is a point in the translated axis system. r is the same point in the original axis system. R = | X   Y   Z | is the position vector of the new origin in the original axis system. Combining eqs (1.1) and (3.1) yields the transformation of the dipole moment under translation of the axis system

 m(1) = ∫ρ(r) (r − R) d3r = ∫ρ(r) r d3r − R ∫ρ(r) d3r = m(1) − R m(0) (3.2)

Similarly, the transformation of the unabridged cartesian moment tensor operators of order 2 through 4 under translation of the axis system is given by

 r(2) = (r − R)⊗(r − R) = r⊗r − R⊗r − r⊗R + R⊗R (3.3) r(3) = r⊗r⊗r − r⊗R⊗r − r⊗r⊗R − R⊗r⊗r + r⊗R⊗R + R⊗R⊗r + R⊗r⊗R − R⊗R⊗R r(4) = r⊗r⊗r⊗r − r⊗r⊗R⊗r − r⊗r⊗r⊗R − r⊗R⊗r⊗r − R⊗r⊗r⊗r + r⊗r⊗R⊗R + r⊗R⊗R⊗r + r⊗R⊗r⊗R + R⊗r⊗R⊗r + R⊗r⊗r⊗R + R⊗R⊗r⊗r − r⊗R⊗R⊗R − R⊗r⊗R⊗R − R⊗R⊗R⊗r − R⊗R⊗r⊗R + R⊗R⊗R⊗R

or, in index notation:

 [r(2)]ij = rirj − riRj − rjRi + RiRj (3.4) [r(3)]ijk = rirjrk − rirjRk − rirkRj − rjrkRi + riRjRk + rjRiRk + rkRiRj − RiRjRk [r(4)]ijkl = rirjrkrl − rirjrkRl − rirjrlRk − rirkrlRj − rjrkrlRi + rirjRkRl + rirlRjRk + rirkRjRl + rjrlRiRk + rjrkRiRl + rkrlRiRj − riRjRkRl − rjRiRkRl − rkRiRjRl − rlRiRjRk + RiRjRkRl ∀ i, j, k, l = x, y, z

Substituting expressions (3.4) for the moment operator in Eq. (1.1) and integrating yields:

 mij = mij − miRj − mjRi + RiRjm0 (3.5) mijk = mijk − mijRk − mikRj − mjkRi + miRjRk + mjRiRk + mkRiRj − RiRjRkm0 mijkl = mijkl − mijkRl − mijlRk − miklRj − mjklRi + mijRkRl + milRjRk + mikRjRl + mjlRiRk + mjkRiRl + mklRiRj − miRjRkRl − mjRiRkRl − mkRiRjRl − mlRiRjRk + RiRjRkRlm0 ∀ i, j, k, l = x, y, z

The results above can be put together into a set of coupled equations, that in matrix notation have the form:

 M = S M (3.6)

where M is the polytensor defined by Eq. (3.1), which contains the original multipole moments, M contains the translated multipole moments and the shift matrix S is a lower unitriangular matrix of dimension N×N, with N = (n+1)(n+2)(n+3)/6, which is factorized as shown in Eq. (3.7):

 m(0) m(1) m(2) m(3) m(4) ⋮
=
 1 01×34 … SD,3×1 I3 03×31 … SQ,6×4 I6 06×25 … SO,10×10 I10 010×15 … SH,15×20 I15 … ⋮ ⋮
 m(0) m(1) m(2) m(3) m(4) ⋮
(3.7)

where the I(n+1)(n+2)/2 are identity matrices, the 0[(n+1)(n+2)/2]×[N−(n+1)(n+2)(n+3)/6] are rectangular matrices whose elements are equal to zero, and the elements of the rectangular matrices S[(n+1)(n+2)/2]×[n(n+1)(n+2)/6] are simple functions of the elements of R:

 SD = − R (3.8)
SQ =
 XX −2X 0 0 XY −Y −X 0 YY 0 −2Y 0 XZ −Z 0 −X YZ 0 −Z −Y ZZ 0 0 −2Z
(3.9)
SO =
 −XXX 3XX 0 0 −3X 0 0 0 0 0 −YYY 0 3YY 0 0 0 −3Y 0 0 0 −ZZZ 0 0 3ZZ 0 0 0 0 0 −3Z −XYY YY 2XY 0 0 −2Y −X 0 0 0 −YXX 2XY XX 0 −Y −2X 0 0 0 0 −ZXX 2XZ 0 XX −Z 0 0 −2X 0 0 −XZZ ZZ 0 2XZ 0 0 0 −2Z 0 −X −YZZ 0 ZZ 2YZ 0 0 0 0 −2Z −Y −ZYY 0 2YZ YY 0 0 −Z 0 −2Y 0 −XYZ YZ XZ XY 0 −Z 0 −Y −X 0
(3.10)
SH =
 XXXX −4XXX 0 0 6XX 0 0 0 0 0 −4X 0 0 0 0 0 0 0 0 0 YYYY 0 −4YYY 0 0 0 6YY 0 0 0 0 −4Y 0 0 0 0 0 0 0 0 ZZZZ 0 0 −4ZZZ 0 0 0 0 0 6ZZ 0 0 −4Z 0 0 0 0 0 0 0 XXXY −3YXX −XXX 0 3XY 3XX 0 0 0 0 −Y 0 0 0 −3X 0 0 0 0 0 XXXZ −3ZXX 0 −XXX 3XZ 0 0 3XX 0 0 −Z 0 0 0 0 −3X 0 0 0 0 YYXY −YYY −3XYY 0 0 3YY 3XY 0 0 0 0 −X 0 −3Y 0 0 0 0 0 0 YYYZ 0 −3ZYY −YYY 0 0 3YZ 0 3YY 0 0 −Z 0 0 0 0 0 0 −3Y 0 ZZXZ −ZZZ 0 −3XZZ 0 0 0 3ZZ 0 3XZ 0 0 −X 0 0 0 −3Z 0 0 0 ZZYZ 0 −ZZZ −3YZZ 0 0 0 0 3ZZ 3YZ 0 0 −Y 0 0 0 0 −3Z 0 0 XXYY −2XYY −2YXX 0 YY 4XY XX 0 0 0 0 0 0 −2X −2Y 0 0 0 0 0 XXZZ −2XZZ 0 −2ZXX ZZ 0 0 4XZ 0 XX 0 0 0 0 0 −2Z −2X 0 0 0 YYZZ 0 −2YZZ −2ZYY 0 0 ZZ 0 4YZ YY 0 0 0 0 0 0 0 −2Y −2Z 0 XXYZ −2XYZ −ZXX −YXX YZ 2XZ 0 2XY XX 0 0 0 0 0 −Z −Y 0 0 0 −2X YYXZ −ZYY −2XYZ −XYY 0 2YZ XZ YY 2XY 0 0 0 0 −Z 0 0 0 0 −X −2Y ZZXY −YZZ −XZZ −2XYZ 0 ZZ 0 2YZ 2XZ XY 0 0 0 0 0 0 −Y −X 0 −2Z
(3.11)

Monopoles, m(0), do not change under translation of the axis system, because they are scalars. Eq. (3.1) shows that also the components of the dipole remain unchanged under translation of the axis system provided that the monopole is null. In general, according to Eqs. (3.5), a symmetric Cartesian tensor m(n) of order n is invariant under translation of the axis system provided that all the symmetric Cartesian tensors of lower order, up to n−1, are null.

On the other hand, the translation of a polytensor, whose component tensors are all null but the monopole m(0), generates a polytensor with non-null component tensors of any order, i.e.
m(1) = −Rm(0),
m(2) = |X2  XY  Y2  XZ  YZ  Z2| m(0),
m(3) = − |X3  Y3  Z3  XY2  YX2  ZX2  XZ2  YZ2  ZY2  XYZ| m(0),   etc.

## 4. - Rotation of Unabridged Moments

In general, a tensor of order n is a mathematical object with n suffixes, Tijk…, which obeys the transformation law

 T'ijk… =   ∑ pqr… LipLjqLkr … Tpqr… (4.1)

where L is the rotation matrix between frames. A scalar, T, is a tensor of order 0, because it is the same in all frames (T' = T). Any vector, like the position vector r, is a tensor of order 1, because r'i = ∑p Lip rp, or in matrix notation r' = L r. For second order tensors, such as the quadrupole tensor, the transformation law T'ij = ∑pq LipLjq Tpq, can be rewritten in matrix notation as T' = L T L, in addition to the tensor notation  T'(2) = L(2,2) T(2).

In the previous sections we defined multidimensional arrays as tensor. Here we define tensors through the transformation low that they obey, as in Eq. (4.1). Definition (4.1) is based on a synecdoche, since the entire tensor is indicated by one of its components (the so called suffix notation). In the present development, the common convention of implied summation over repeated indices is not being used. If we did, we would not bother to write down the ∑ in the right side of Eq. (4.1) because suffixes p, q and r appear twice in a single term of an expression.

As an example, let's write Eq. (4.1) for the second order Cartesian moment m(2), using the suffix notation:

 m'ij = ∑ q=x,y,z ∑ p=x,y,z LipLjqmpq = = LixLjxmxx + LiyLjxmyx + LizLjxmzx + LixLjymxy + LiyLjymyy + LizLjymzy + LixLjzmxz + LiyLjzmyz + LizLjzmzz (4.2)

where the two sums in the second member have been developed following the algorithm of Program 1, so that the nine tensor components appear in canonical order in the sum of the last member. Since Cartesian tensors, m(n), are totally symmetric, there are only (n+1)(n+2)/2 distinct components (namely, six for n = 2). The number of times that each of the distinct components appears in the complete set of 3n components is given by Eq. (2.2). Then, Eq. (4.2) can be rewritten as

 m'l(2) = (n+1)(n+2)/2 ∑ k=1 mk(2) ( gk ∑ {pq} LipLjq ) = = LixLjxmxx + (LiyLjx + LixLjy) myx + (LizLjx + LixLjz) mzx + + LiyLjymyy + (LizLjy + LiyLjz) mzy + LizLjzmzz (4.3)

where k and l span the compressed canonical array which comprises (n+1)(n+2)/2 components; {pq} spans the gk doublets (p,q) which have mk(2) as a common factor, and (i,j) is the doublet of suffixes which corresponds to the l-th suffix of the compressed canonical array.

The last member of Eq. (4.2) defines, by elements, the rotation matrix L(2,2) of the complete canonical tensor m(2):

L(2,2) = LL =
 LxxL LxyL LxzL LyxL LyyL LyzL LzxL LzyL LzzL
=

 LxxLxx LxyLxx LxzLxx LxxLxy LxyLxy LxzLxy LxxLxz LxyLxz LxzLxz LyxLxx LyyLxx LyzLxx LyxLxy LyyLxy LyzLxy LyxLxz LyyLxz LyzLxz LzxLxx LzyLxx LzzLxx LzxLxy LzyLxy LzzLxy LzxLxz LzyLxz LzzLxz LxxLyx LxyLyx LxzLyx LxxLyy LxyLyy LxzLyy LxxLyz LxyLyz LxzLyz LyxLyx LyyLyx LyzLyx LyxLyy LyyLyy LyzLyy LyxLyz LyyLyz LyzLyz LzxLyx LzyLyx LzzLyx LzxLyy LzyLyy LzzLyy LzxLyz LzyLyz LzzLyz LxxLzx LxyLzx LxzLzx LxxLzy LxyLzy LxzLzy LxxLzz LxyLzz LxzLzz LyxLzx LyyLzx LyzLzx LyxLzy LyyLzy LyzLzy LyxLzz LyyLzz LyzLzz LzxLzx LzyLzx LzzLzx LzxLzy LzyLzy LzzLzy LzxLzz LzyLzz LzzLzz
(4.4)

On the other hand, the rotation matrix of the corresponding compressed canonical tensor is defined by the last member of Eq. (4.3):

L(2,2) = C L(2,2) C g(2) =

 Lxx2 2 LxxLxy 2 LxxLxz Lxy2 2 LxyLxz Lxz2 LyxLxx LyyLxx + LyxLxy LyzLxx + LyxLxz LyyLxy LyzLxy + LyyLxz LyzLxz LzxLxx LzyLxx + LzxLxy LzzLxx + LzxLxz LzyLxy LzzLxy + LzyLxz LzzLxz Lyx2 2 LyxLyy 2 LyxLyz Lyy2 2 LyyLyz Lyz2 LzxLyx LzyLyx + LzxLyy LzzLyx + LzxLyz LzyLyy LzzLyy + LzyLyz LzzLyz Lzx2 2 LzxLzy 2 LzxLzz Lzy2 2 LzyLzz Lzz2
(4.5)

where the diagonal matrix g(n) and the rectangular matrix C are defined by Eqs (2.2-5).

It's worth noting that, given the rotated position vector r' = L r, the primed moment operator of order 2, r'(2), according to Eq. (1.4), is given by the Kronecker product r'r', which leads to:

 r'(2) = r'⊗r' = (L r)⊗(L r) = (L⊗L)(r⊗r) = L(2,2) r(2) (4.6)

where the mixed-product property has been applied.[2Zhang, H.; Ding, F.
Journal of Applied Mathematics 2013, Article ID 296185, 8 pages.
] Eq. (4.6) can be generalized to r'(n) = L(n,n)r(n) and integrated by Eq. (1.1), yielding

 m'(n) = L(n,n) m(n) (4.7)

## 5. - Tensor Contraction and Traces

The contraction of a tensor is obtained by setting unlike indices equal and summing. Contraction reduces the tensor order by 2. Contraction is not defined for tensors of order 0 (scalars) and order 1 (vectors). For a second order tensor m(2), which has two indices, the contraction is the same as the trace (provided that m(2) is interpreted as a matrix) and is therefore a scalar:

 m(2:1) = ∑ i=x,y,z mii = mxx + myy + mzz ≡ Tr [m(2)] (5.1)

For a totally symmetric third order cartesian tensor m(3), which has three indices, contraction implies three trace relationships, each one defining a component of a first order tensor or vector, m(3:1):

m(3:1) =
 mx(3:1) my(3:1) mz(3:1)
 m100(3:1) m010(3:1) m001(3:1)
 ∑ i=x,y,z miix ∑ i=x,y,z miiy ∑ i=x,y,z miiz
=
 mxxx + myyx + mzzx myyy + mxxy + mzzy mzzz + mxxz + myyz
(5.2)

Finally, for a totally symmetric fourth order cartesian tensor m(4), which has four indices, contraction implies six distinct trace relationships, each one defining a component of a symmetric compressed second order tensor, m(4:1):

m(4:1) =
 mxx(4:1) mxy(4:1) myy(4:1) mxz(4:1) myz(4:1) mzz(4:1)
 m200(4:1) m110(4:1) m020(4:1) m101(4:1) m011(4:1) m002(4:1)
 ∑ i=x,y,z miixx ∑ i=x,y,z miixy ∑ i=x,y,z miiyy ∑ i=x,y,z miixz ∑ i=x,y,z miiyz ∑ i=x,y,z miizz
=
 mxxxx + myyxx + mzzxx mzzxy + mxxxy + myyxy myyyy + mxxyy + mzzyy myyxz + mxxxz + mzzxz mxxyz + myyyz + mzzyz mzzzz + mxxzz + myyzz
(5.3)

The double contraction of a tensor is obtained by setting four unlike indices two by two to be equal and summing over the two pairs of equal indices. Double contraction reduces the tensor order by 4. Therefore, double contraction of a totally symmetric fourth order cartesian tensor m(4) yields a scalar:

 m(4:2) = m000(4:2) = ∑ i=x,y,z ∑ j=x,y,z miijj = ∑ i=x,y,z mii(4:1) = Tr [ m(4:1) ] = mxxxx + myyyy + mzzzz + 2 (mxxyy + mxxzz + myyzz) (5.4)

In general,[12Applequist, Jon
Chem. Phys. 1984, 85, 279-290.
; 13Applequist, Jon
J. Chem. Phys. 1985, 83, 809-826.
; 14Applequist, Jon
J. Phys. A: Math. Gen. 1989, 22, 4303-4330.
] the trace of a tensor m(n) of order n with respect to one pair of suffixes is a tensor m(n:1) of order n−2 and is given by elements:

 m(n:1)r3r4…rn = ∑ i=x,y,z m(n)iir3r4…rn = m(n)kx+2,ky,kz + m(n)kx,ky+2,kz + m(n)kx,ky,kz+2 = m(n:1)kx,ky,kz (5.5)

where kx + ky + kz = n − 2. The distinct components of m(n:1) are exactly n(n−1)/2, one for each selection of the (n−2) suffixes x, y and z. The trace of a tensor m(n) of order n with respect to l pairs of suffixes (called an “l-fold trace”) is a tensor m(n:l) of order (n − 2l) and is given by elements:[14Applequist, Jon
J. Phys. A: Math. Gen. 1989, 22, 4303-4330.
]

 m(n:l)r2l+1…rn = ∑ r1=x,y,z   ∑ r2=x,y,z … ∑ rl=x,y,z m(n)r1r1 r2r2 … rlrl r2l+1 … rn = ∑ lxlylz l! / lx! ly! lz! m(n)kx+2lx,ky+2ly,kz+2lz = m(n:l)kx,ky,kz (5.6)

where kx + ky + kz = n − 2l and the sum is over all non-negative indices such that lx + ly + lz = l. The trinomial coefficient g(l;lxlylz) = l!/lx!ly!lz! appears here as the number of the ways one can place lx pairs xx, ly pairs yy, and lz pairs zz in the indices r1r1 r2r2rlrl of the l-fold trace.[14Applequist, Jon
J. Phys. A: Math. Gen. 1989, 22, 4303-4330.
] For instance, if l = 2 there are six distinct sets of indices lxlylz (namely, 200, 020, 002, 110, 101, 011, with multiplicity 1, 1, 1, 2, 2, 2, respectively) so that expression (5.4) is easily recovered. If l = 1 there are three distinct sets of indices, i.e. 100, 010, 001, so that eqs (5.4), (5.3) and (5.2) are obtained.

If the trace vanishes regardless of which index pair is contracted, the tensor is said to be totally traceless. If a tensor is totally symmetric and traceless in one index pair, then it is traceless for all index pairs, and is said to be totally symmetric and traceless.

## 6. - Traceless Moments

The traceless cartesian moment tensor of order n of a charge distribution ρ(r) is defined by elements in component index notation[15Stone, A. J.
"The Theory of Intermolecular Forces"
Oxford University Press, Oxford, 2000; p. 16.
; 16Coppens, P.
"X-Ray Charge Densities and Chemical Bonding"
International Union of Crystallography, Oxford University Press, Oxford, 1997; p. 144.
]

 μ r1r2…rn−1rn = (−1)n / n! ∫ρ(r) r2n+1 ∂n r−1 / ∂r1∂r2…∂rn−1∂rn d3r ∀ ri = x, y, z   ∧   i = 1, 2, …, n (7.1)

following the convenction that ri subscripts denote cartesian axes x, y, z, or by elements in degree index notation[14Applequist, Jon
J. Phys. A: Math. Gen. 1989, 22, 4303-4330.
]

 μnxnynz = (−1)n / n! ∫ρ(r) r2n+1 ∂nx+ny+nz (x2 + y2 + z2)−½ / ∂xnx ∂yny ∂znz d3r ∀ nx, ny, nz = 0, 1, 2, …, n   ∧   nx + ny + nz = n (7.2)

where ni is the number of times i occurs in the index set r1r2rn−1rn and nx + ny + nz = n. The latter notation is usefull for compressed tensors, as it addresses only distinct components of the complete tensors.

Direct differentiation of r−1 yields in index notation

 ∂ r−1 / ∂ ri = − r−3 ri (7.3) ∂2 r−1 / ∂ ri ∂ rj = r−5 [ 3 rirj − r2 δij ] ∂3 r−1 / ∂ ri ∂ rj ∂ rk = − r−7 [ 15 ri rj rk − 3 r2 (ri δjk + rj δik + rk δij )] ∂4 r−1 / ∂ ri ∂ rj ∂ rk ∂ rl = r−9 [ 105 ri rj rk rl − 15 r2 (ri rj δkl + ri rk δjl + ri rl δjk + rj rk δil + rj rl δik + rk rl δij) + 3 r4 (δijδkl + δikδjl + δilδjk) ] … ∂n r−1 / ∂r1 ∂r2 … ∂rn−1 ∂rn = (−1)n r−(2n+1) ⌊n/2⌋ ∑ m=0 (−1)m(2n − 2m − 1)!! ×   ∑ T{ri} δr1r2 δr3r4 … δr2m−1r2m [r(n:m)] r2m+1 … rn

where ⌊n/2⌋ denotes the integer part of n/2; (2n − 2m − 1)!! = (2n − 2m)!/2n−m (n − m)! = 1 ⋅ 3 ⋅ 5 ⋅ (2n − 2m − 1) with (−1)!! = 1; δij is the Kronecker delta; and the sum over T{ri} is the sum over all permutations of the symbols r1 r2rn which give distinct terms.

Inserting the moment operators of Eq. (7.3) into Eq. (7.1) and integrating yields:[12Applequist, Jon
Chem. Phys. 1984, 85, 279-290.
]

 μi = mi (7.9) Θij = 3⁄2 mij − 1⁄2 ⟨r2⟩ δij (7.10) Ωijk = 5⁄2 mijk − 1⁄2 ⟨r2⟩ (mi δjk + mj δik + mk δij) (7.11) Φijkl = 35⁄8 mijkl − 5⁄8 ⟨r2⟩ (mij δkl + mik δjl + mil δjk + mjk δil + mjl δik + mkl δij) + 1⁄8 ⟨r4⟩ (δij δkl + δik δjl + δil δjk) (7.12)

The order n of tensors m has been omitted in Eqs (7.9-12) because its value is clear from the number of suffixes. The symbols μ, Θ, Ω, and Φ, have been adopted [17Buckingham, A. D.
J. Chem. Phys. 1959, 30, 1580-1585.
; 18Buckingham, A. D.
Advan. Chem. Phys. 1967, 12, 107-142.
; 19McLean, A. D.; Yoshimine, M. J. Chem. Phys. 1967, 47, 1927-1935.] for μ(1), μ(2), μ(3), and μ(4), respectively. In addition:

 ⟨r2⟩ = ⟨x2 + y2 + z2⟩ = mxx + myy + mzz (7.13) ⟨r4⟩ = ⟨x4 + y4 + z4 + 2 (x2y2 + x2z2 + y2z2)⟩ = mxxxx + myyyy + mzzzz + 2 (mxxyy + mxxzz + myyzz) (7.14)

According to Eq. (5.6), the elements of tensor r(n:m), which is the trace of r(n) with respect to m pairs of suffixes (called an “m-fold trace”), can be written explicitly as

 [r(n:m)]r2m+1…rn = ∑ r1r2…rm [r(n)] r1r1 r2r2 … rmrm r2m+1 … rn= ∑ r1=x,y,z r12 ∑ r2=x,y,z r22 … ∑ rm=x,y,z rm2 r2m+1 … rn =   r2m r2m+1 … rn =   r2m xkx yky zkz (7.4)

where kx + ky + kz = n − 2m.   In addition[21Trinomial expansion and trinomial coefficients]

 r2m   =   (x2 + y2 + z2)m   = ∑ mx+my+mz = m m! / mx! my! mz! x2mx y2my z2mz (7.4)

so that

 [r(n:m)]r2m+1…rn = ∑ mx+my+mz = m m! / mx! my! mz! xkx+2mx yky+2my zkz+2mz (7.4)

where the sum is over all non-negative indices mx, my, mz such that mx + my + mz = m,  and 2miki  ∀ i = x,y,z,  and kx + ky + kz = n − 2m.   [Please, note that in the equations above ri denotes any cartesian coordinate x, y, z, whereas r without suffix denotes |r| = (x2 + y2 + z2)−½].

In tensor notation, equations (7.1) and (7.2) become:

 μ(n) = (−1)n / n! ∫ρ(r) r2n+1 ∇ (n) r−1 d3r (7.1)

where the nabla symbol denotes the vector differential operator:

= / r =
 ∂ / ∂ x ∂ / ∂ y ∂ / ∂ z
(7.2)

and 2 denotes the Hessian matrix operator

2 = = vec ( ) = vec
 ∂2 / ∂ x2 ∂2 / ∂ x ∂ y ∂2 / ∂ x ∂ z ∂2 / ∂ y ∂ x ∂2 / ∂ y2 ∂2 / ∂ y ∂ z ∂2 / ∂ z ∂ x ∂2 / ∂ z ∂ y ∂2 / ∂ z2
(7.2)

whose trace is the laplacian operator, 2:

 ∇(2:1) = Tr (∇2) = ∇ ⋅ ∇ = ∇† ∇2 = ∇ 2 = ∂2 / ∂ x2 + ∂2 / ∂ y2 + ∂2 / ∂ z2 (7.2)

The traceless condition with respect to any pair of suffixes comes from the fact that the laplacian of r−1 is null:[15Stone, A. J.
"The Theory of Intermolecular Forces"
Oxford University Press, Oxford, 2000; p. 16.
]

 ∇ 2 r−1 = 0. (7.2)

The last equation of system (7.3) shows that the n-th gradient of r−1 is a linear combination of r(n) and its traces in all pairs of indices.[12Applequist, Jon
Chem. Phys. 1984, 85, 279-290.
] In other words, the totally symmetric tensor r(n) of order n can be transformed in a traceless totally symmetric tensor (n)r−1 by the action of the detracer operator 𝔇n,[12Applequist, Jon
Chem. Phys. 1984, 85, 279-290.
; 13Applequist, Jon
J. Chem. Phys. 1985, 83, 809-826.
] which linearly combines the tensor components through the associated square matrix D(n,n):[14Applequist, Jon
J. Phys. A: Math. Gen. 1989, 22, 4303-4330.
]

 ∇(n)r−1 = (−1)n r−(2n+1) D(n,n)r(n) (7.5)

or in index notation

 [∇(n)r−1] β1 β2 … βn = (−1)n r−(2n+1) ∑ α1α2…αn [D(n,n)] β1 β2 … βn, α1 α2 … αn [r(n)] α1 α2 … αn (7.6)

Inserting Eq. (7.5) into Eq. (7.1) and integrating, yields:

 μ(n) = 1 / n! D(n,n)m(n) (7.7)

From (7.3) and (7.4) it can be seen that the components of the tensor of subdivided order D(n,n) [3Applequist, J.
J. Math. Phys. 1983, 24, 736-741.
] have the following expression:[14Applequist, Jon
J. Phys. A: Math. Gen. 1989, 22, 4303-4330.
]

 [D(n,n)] β1 β2 … βn, α1 α2 … αn = ⌊n/2⌋ ∑ m=0 (−1)m(2n − 2m − 1)!! ×   ∑ T{αβ} δα1α2 … δα2m−1α2m δβ1β2 … δβ2m−1β2m δα2m+1β2m+1 … δαnβn (7.8)

where the sum over T{αβ} is the sum over all permutations of the symbols α1 α2αn and β1 β2βn giving distinct terms. For example, Eq. (7.8) gives for n = 2

 [D(2,2)]ij,kl = δik δjl − ⅓ δij δkl
D(2,2) =
 2 0 −1 0 0 −1 0 3 0 0 0 0 −1 0 2 0 0 −1 0 0 0 3 0 0 0 0 0 0 3 0 −1 0 −1 0 0 2
(7.2)
 [D(3,3)] β1 β2 β3, α1 α2 α3 = 15 δα1β1 δα2β2 δα3β3 − 3 ( δα1α2 δβ1β2 δα3β3 + δα1α3 δβ1β3 δα2β2 + δα2α3 δβ2β3 δα1β1 )

which, using Eq. (7.6), gives the third equation of system (7.3).

Table − Using the last equation of system (7.3) to find the first four equations of the same system.
α1 α2αn n m (−1)m(2n − 2m − 1)!!   T{α}
α 1 0 −1 rα
α1α2 2 0 3 rα1 rα2
1 −1 δα1α2
α1α2α3 3 0 15 rα1 rα2 rα3
1 −3 δα1α2 rα3 + δα1α3 rα2 + δα2α3 rα1
α1α2α3α4 4 0 105 rα1 rα2 rα3 rα4
1 −15 δα1α2 rα3rα4 + δα1α3 rα2rα4 + δα1α4 rα2rα3 + δα2α3 rα1rα4 + δα2α4 rα1rα3 + δα3α4 rα1rα2
2 3 δα1α2 δα3α4 + δα1α3 δα2α4 + δα1α4 δα2α3
D(3,3) =
 2 0 0 −3 0 0 −3 0 0 0 0 2 0 0 −3 0 0 −3 0 0 0 0 2 0 0 −3 0 0 −3 0 −1 0 0 4 0 0 −1 0 0 0 0 −1 0 0 4 0 0 −1 0 0 0 0 −1 0 0 4 0 0 −1 0 −1 0 0 −1 0 0 4 0 0 0 0 −1 0 0 −1 0 0 4 0 0 0 0 −1 0 0 −1 0 0 4 0 0 0 0 0 0 0 0 0 0 5
(7.2)
D(4,4) =
 8 3 3 0 0 0 0 0 0 −24 −24 6 0 0 0 3 8 3 0 0 0 0 0 0 −24 6 −24 0 0 0 3 3 8 0 0 0 0 0 0 6 −24 −24 0 0 0 0 0 0 20 0 −15 0 0 0 0 0 0 0 0 −15 0 0 0 0 20 0 0 −15 0 0 0 0 0 −15 0 0 0 0 −15 0 20 0 0 0 0 0 0 0 0 −15 0 0 0 0 0 0 20 0 −15 0 0 0 −15 0 0 0 0 0 0 −15 0 0 20 0 0 0 0 0 −15 0 0 0 0 0 0 0 −15 0 20 0 0 0 −15 0 0 −4 −4 1 0 0 0 0 0 0 27 −3 −3 0 0 0 −4 1 −4 0 0 0 0 0 0 −3 27 −3 0 0 0 1 −4 −4 0 0 0 0 0 0 −3 −3 27 0 0 0 0 0 0 0 0 0 −5 0 −5 0 0 0 30 0 0 0 0 0 0 −3 0 0 −3 0 0 0 0 0 30 0 0 0 0 −3 0 −3 0 0 0 0 0 0 0 0 30
(7.2)

## References and Notes

1. If A is an m × n matrix and B is a p×q matrix, then the Kronecker product AB is the mp×nq block matrix:
AB
 A11B … A1nB ⋮ ⋱ ⋮ Am1B ⋯ AmnB

See reference 2 and external links L7 and L8 for more details.
2. "On the Kronecker Products and Their Applications"
Zhang, H.; Ding, F.   Journal of Applied Mathematics 2013, Article ID 296185, 8 pages,
http://dx.doi.org/10.1155/2013/296185
3. "Cartesian polytensors"
Applequist, J. J. Math. Phys. 1983, 24, 736-741.
4. 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 column-major 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.
5. A dyadic product is the special case of the Kronecker product or tensor product between two vectors of the same dimension. If a = |ax ay az| and b = |bx by bz| are two cartesian vectors, their dyadic products are:
ab =
 ax b ay b az b
=
 ax bx ax by ax bz ay bx ay by ay bz az bx az by az bz
=
 bx ax by ax bz ax bx ay by ay bz ay bx az by az bz az
= vec
 bx ax bx ay bx az by ax by ay by az bz ax bz ay bz az
= vec
 bx a† by a† bz a†
= vec (b a)
and, similarly:
 b⊗a = vec (a b†)
The two equations above show that the dyadic product is not commutative. The tensor product of two vectors of different dimensions is called outer product [cf Link 6]. It contrasts with the inner product of vectors a and b, which yields a scalar:
 ⟨a, b⟩ = ⟨b, a⟩ = a†b = b†a = axbx + ayby + azbz
and is better known as the dot product. The inner product is the trace of the outer product.
6. "The commutation matrix: some properties and applications"
Magnus, J. R.; Neudecker, H. The Annals of Statistics 1979, 7, 381-394.
7. "The Vec-Permutation Matrix, The Vec Operator and Kronecker Products: A Review"
Henderson, H. V.; Searle, S. R. Linear and Multilinear Algebra 1981, 9, 271-288.
8. In general, a tensor r(2,0) with components in anticanonical order can be transformed into an equivalent tensor with components in canonical order by means of a suitable vec-permutation matrix,[7Henderson, H. V.; Searle, S. R.
Linear and Multilinear Algebra 1981, 9, 271-288.
] also known as commutation matrix:[6Magnus, J. R.; Neudecker, H.
The Annals of Statistics 1979, 7, 381-394.
]
 xx yx zx xy yy zy xz yz zz
=
 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1
 xx xy xz yx yy yz zx zy zz

9. "Nested summation symbols and perturbation theory"
Carbó, R.; Besalú, E.   J. Math. Chem. 1993, 13, 331-342.
10. "Definition and quantum chemical applications of nested summation symbols and logical functions: Pedagogical artificial intelligence devices for formulae writing, sequential programming and automatic parallel implementation"
Carbó, R.; Besalú, E.   J. Math. Chem. 1995, 18, 37-72.
11. Gaussian 09, Revision E.01, Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian, Inc., Wallingford CT, 2009.
12. "Fundamental relationships in the theory of electric multipole moments and multipole polarizabilities in static fields"
Applequist, Jon Chem. Phys. 1984, 85, 279-290.
13. "A multipole interaction theory of electric polarization of atomic and molecular assemblies"
Applequist, Jon J. Chem. Phys. 1985, 83, 809-826.
14. "Traceless Cartesian tensor forms for spherical harmonic functions: new theorems and applications to electrostatics of dielectric media"
Applequist, Jon J. Phys. A: Math. Gen. 1989, 22, 4303-4330.
15. Stone, A. J. "The Theory of Intermolecular Forces", Oxford University Press, Oxford, 2000, ISBN 0-19-855883-X; p. 16.
16. Coppens, P. "X-Ray Charge Densities and Chemical Bonding", International Union of Crystallography, Oxford University Press, Oxford, 1997, ISBN: 9780195098235; p. 144.
17. "Direct Method of Measuring Molecular Quadrupole Moments"
Buckingham, A. D. J. Chem. Phys. 1959, 30, 1580-1585.
18. "Permanent and Induced Molecular Moments and Long-Range Intermolecular Forces"
Buckingham, A. D. Advan. Chem. Phys. 1967, 12, 107-142.
19. "Theory of Molecular Polarizabilities"
McLean, A. D.; Yoshimine, M. J. Chem. Phys. 1967, 47, 1927-1935.
20. The product C C is not positive definite and cannot be inverted. In addition, C g(n) CI.
21. The n-th power of a trinomial a + b + c can be expanded into a sum of (n+1)(n+2)/2 monomials according to the following equation
 (a + b + c)n = ∑ i+j+k = n ( n   i,j,k ) ai bj ck
where i, j, k are non-negative integers and
 ( n   i,j,k ) = n! / i! j! k!
are the trinomial coefficients.   Their sum over all non-negative indices i, j, k such that i + j + k = n equals 3n:
 ∑ i+j+k = n ( n   i,j,k ) = 3n.
A trinomial coefficient can be expressed as the product of three binomial coefficients:
 ( n   i,j,k ) = ( n   i ) ( n−i   j ) ( n−i−j   k )
This identity can be verified by applying the definition of binomial coefficient, i.e. ( n   i ) = n! / i! (n−i)!. It's worth noting that the binomial coefficient ( n   i ) can also be written as ( n   i,j ) = n! / i! j!, provided that j = n − i.