Houxian Chen, Menglin Liu, Tianying Yan^,Institute of New Energy Material Chemistry, School of Materials Science and Engineering, National Institute for Advanced Materials, Nankai University, Tianjin 300350, China

Received:2020-01-21Revised:2020-03-04Accepted:2020-03-19Online:2020-06-18


Abstract
The Buckingham expansion is important for understanding molecular multipoles and (hyper)polarizabilities. In this study, we give a complete derivation of the Buckingham expansion in the traced form using successive Taylor series. Based on the derivation results, a general Buckingham expansion in the traced form is proposed, from which highly accurate numerical calculations using the finite field method can be achieved. The transformations from the traced multipoles and multipole–multipole polarizabilities to the corresponding traceless counterparts are realized with an auxiliary traced electric field gradient. The applications of the finite field method in this study show good agreements with previous theoretical calculations and experimental measurements.
Keywords： Buckingham expansion;multipole;(hyper)polarizability;finite field method;ab initio calculations


PDF (647KB)MetadataMetricsRelated articlesExportEndNote|Ris|BibtexFavorite
Cite this article
Houxian Chen, Menglin Liu, Tianying Yan. Molecular multipoles and (hyper)polarizabilities from the Buckingham expansion: revisited. Communications in Theoretical Physics, 2020, 72(7): 075503- doi:10.1088/1572-9494/ab8a0d

1. Introduction

Multipoles and (hyper)polarizabilities are important molecular interactions [1, 2]. The interactions between molecules and external electric fields and field gradients can be expanded with multipoles in external electric fields and field gradients, which is the well-known multipole expansion. Furthermore, if the molecule is polarizable, the multipoles are also inducible by the field and field gradients. The resultant interactions are expanded to include all the permanent and induced multipoles, with the induced components manifested in the (hyper)polarizabilities, and we denote in this study the resultant expression as the Buckingham expansion [3–5]. Thus, as multiple external electric fields and field gradients E are applied to the polarizable charge distribution ρ(r), which can be expanded by permanent and induced multipoles M(r), the responses via the Buckingham expansion are as shown schematically in figure 1.

Figure 1.

New window|Download| PPT slide
Figure 1.A schematic Buckingham expansion, in which E can be external electric fields and field gradients, ρ(r) represents the polarizable charge distribution, andM(r) denotes permanent and induced multipoles. The equation corresponds to the Buckingham expansion in the traceless form, which is derived in equation (27).


The Buckingham expansion is not only of interest theoretically, but also leads to important applications in determining molecular (hyper)polarizabilities [6] and quadrupoles [3, 7]. By designing a four-wire electric field gradient condenser, which induces birefringence (Kerr effects) on the polar molecular system, inherently originated from the anisotropic molecular multipoles/(hyper)polarizabilities as represented by the Buckingham expansion, Buckingham derived the method to detect the molecular quadrupole [3]. This theoretical derivation was realized experimentally with direct detection of the quadrupole of carbon dioxide [7].

Since it is generally difficult to detect molecular hyperpolarizabilities experimentally, except in some small molecules of high symmetry, it is common nowadays to calculate molecular multipoles and (hyper)polarizabilities numerically against high-level ab initio calculations, either using test charges to generate the desired finite external electric field and field gradients [8, 9], or manipulating the electric field tensors directly [10, 11]. Either way, the molecular multipoles and (hyper)polarizabilities can be obtained by finite difference with high accuracy up to O(Eⁿ), in which E denotes the external fields or field gradients. In order to get numerical multipoles and (hyper)polarizabilities, especially for the multipole–multipole polarizabilities, accurately up to high-order n, the complete Buckingham expansion needs to be known.

To the best of our knowledge, the Buckingham expansion was originally proposed up to the fourth order of electric fields or field gradients in the traceless form, in which all the multipoles and (hyper)polarizabilities are traceless [3, 4]. For ab initio calculations, it is more convenient to work instead with traced multipoles and (hyper)polarizabilities. Thus, a general Buckingham expansion in the traced form up to arbitrary order of electric fields or field gradients, as well as the corresponding transformation from traced to traceless multipoles and (hyper)polarizabilities, are desired. Aiming to fill in the gaps, McLean and Yoshimine [12], Applequist [13], and Pedersen et al [11] all derived the transformation from traced to traceless multipoles and (hyper)polarizabilities, in which the transformation of the multipole–multipole polarizabilities were achieved by comparisons between the Buckingham expansion in the traced and traceless forms. Considering the importance of the Buckingham expansion, we believe that a revisit of it is needed in order to clarify some of its ambiguities. In this study, we give a complete derivation of the Buckingham expansion, which we think is easily understood and more closely resembles its original form, as we revisit this important and interesting topic.

In this study, the derivation is based on solely Cartesian tensors, because they are more relevant to the experimentally reported molecular multipoles and (hyper)polarizabilities. The transformation between Cartesian representation and spherical representation with solid harmonic spherical bases has been well documented [2, 11]. In section 2 we derive the Buckingham expansion in the traced form, using Taylor series, with necessary discussions on the Buckingham convention according to the original Buckingham expansion in the traceless form. We next propose an auxiliary traced field gradient, from which we derive the Buckingham expansion in the traceless form in section 3. By using this method, we can easily obtain the transformation from traced multipoles and (hyper)polarizabilities to their corresponding traceless counterparts. While the former is well known, the latter may not be trivial without the auxiliary traced field gradient. Additionally, we demonstrate that the traceless features of multipoles as well as the multipole–multipole polarizabilities, are consequences of the Laplace equation. In section 4, we apply the finite field method to high-level ab initio calculations, utilizing the Buckingham expansion, to numerically calculate the dipole polarizability, quadrupole and quadrupole–quadrupole polarizability of carbon dioxide. In appendix A, we briefly summarize previous derivations, especially McLean and Yoshimine’s derivation [12], and derive the conversions of the multipole–multipole polarizabilities of different conventions. In appendix B, we propose a general Buckingham expansion in the traced form and derive the finite difference expressions of the traced quadrupole and quadrupole–quadrupole polarizability with accuracy of O(E⁴). Note that the derivation of any multipole or (hyper)polarizability is straightforward using the finite field method in appendix B, and the transformations from the traced quadrupole and quadrupole–quadrupole polarizability to the corresponding traceless counterparts are derived in section 3. In appendix C, we propose a general formula for the Buckingham expansion in the traceless form based on the discussions in the text.

2. The Buckingham expansion in the traced form

Considering the interaction of a charge distribution ρ(r), centered at the origin O of a Cartesian coordinate system, with an external electrostatic perturbation, the external electrostatic potential at r can be expanded via Taylor series around O, i.e. [2, 12]:(1)$ \begin{eqnarray}\begin{array}{lll}\varphi ({\boldsymbol{r}}) & = & \varphi (0)-{r}_{\alpha }{E}_{\alpha }-\displaystyle \frac{1}{2}{r}_{\alpha }{r}_{\beta }{E}_{\alpha \beta }-\displaystyle \frac{1}{6}{r}_{\alpha }{r}_{\beta }{r}_{\gamma }{E}_{\alpha \beta \gamma }\\ & & -\,\displaystyle \frac{1}{24}{r}_{\alpha }{r}_{\beta }{r}_{\gamma }{r}_{\delta }{E}_{\alpha \beta \gamma \delta }-\cdots ,\end{array}\end{eqnarray}$in which φ, E_α, E_αβ, E_αβγ and E_αβγδ denote external electric potential, field, gradients, etc, and the Greek subscripts denote the Cartesian coordinate component $\{x,y,z\}.$ In equation (1), Einstein summation convention on the repeated subscripts is adopted throughout.

Using equation (1), the interaction between ρ(r) and the externally perturbing φ(r) is given by [2]:(2)$ \begin{eqnarray}\begin{array}{lll}U-{U}_{0} & = & \displaystyle {\int }_{V}{\rm{d}}v\rho ({\boldsymbol{r}})\varphi ({\boldsymbol{r}})=q\varphi -{\mu }_{\alpha }^{{\prime} }{E}_{\alpha }\\ & & -\,\displaystyle \frac{1}{2}{Q}_{\alpha \beta }^{{\prime} }{E}_{\alpha \beta }-\displaystyle \frac{1}{6}{O}_{\alpha \beta \gamma }^{{\prime} }{E}_{\alpha \beta \gamma }\\ & & -\,\displaystyle \frac{1}{24}{H}_{\alpha \beta \gamma \delta }^{{\prime} }{E}_{\alpha \beta \gamma \delta }-\cdots ,\end{array}\end{eqnarray}$in which U₀ is the energy without external field perturbations, and q, ${\mu }_{\alpha }^{{\prime} },$${Q}_{\alpha \beta }^{{\prime} },$${O}_{\alpha \beta \gamma }^{{\prime} },$ and ${H}_{\alpha \beta \gamma \delta }^{{\prime} }$ are monopole (charge), dipole, traced quadrupole, traced octopole and traced hexadecapole, respectively, i.e.:(3a)$ \begin{eqnarray}q=\displaystyle {\int }_{V}{\rm{d}}v\rho ({\boldsymbol{r}})\end{eqnarray}$(3b)$ \begin{eqnarray}{\mu }_{\alpha }^{{\prime} }=\displaystyle {\int }_{V}{\rm{d}}v\rho ({\boldsymbol{r}}){r}_{\alpha }\end{eqnarray}$(3c)$ \begin{eqnarray}{Q}_{\alpha \beta }^{{\prime} }=\displaystyle {\int }_{V}{\rm{d}}v\rho ({\boldsymbol{r}}){r}_{\alpha }{r}_{\beta }\end{eqnarray}$(3d)$ \begin{eqnarray}{O}_{\alpha \beta \gamma }^{{\prime} }=\displaystyle {\int }_{V}{\rm{d}}v\rho ({\boldsymbol{r}}){r}_{\alpha }{r}_{\beta }{r}_{\gamma }\end{eqnarray}$(3e)$ \begin{eqnarray}{H}_{\alpha \beta \gamma \delta }^{{\prime} }=\displaystyle {\int }_{V}{\rm{d}}v\rho ({\boldsymbol{r}}){r}_{\alpha }{r}_{\beta }{r}_{\gamma }{r}_{\delta }\end{eqnarray}$
$ \begin{eqnarray*}\vdots \end{eqnarray*}$Equation (3) stands because the external electrostatic perturbations, φ, E_α, E_αβ, E_αβγ, E_αβγδ, etc, are independent of r and the integration is performed over the entire space V occupied by ρ(r), which is well separated from the external source. The primes on the superscripts of the multipoles in equations (2) and (3) remind us that they may all be inducible if ρ(r) is polarizable, i.e.:(4)$ \begin{eqnarray}\begin{array}{l}{\mu }_{\alpha }^{{\prime} }={\mu }_{\alpha }+{\mu }_{\alpha }^{{\rm{ind}}},{Q}_{\alpha \beta }^{{\prime} }={Q}_{\alpha \beta }+{Q}_{\alpha \beta }^{{\rm{ind}}},\\ {O}_{\alpha \beta \gamma }^{{\prime} }={O}_{\alpha \beta \gamma }+{O}_{\alpha \beta \gamma }^{{\rm{ind}}},\,\cdots ,\end{array}\end{eqnarray}$in which μ_α, Q_αβ, and O_αβγ are the permanent multipoles, with the superscript ‘ind’ in the second term on the RHS denotes the individually inducible counterpart. Such classification is subtle, because the inducible terms only appear with additional external perturbations. Otherwise, they are the permanent multipoles. Since there is no charge source or sink other than ρ(r), q is permanent and not inducible.

With additional external perturbation, the multipoles can be expanded, again, via Taylor series, i.e.:(5a)$ \begin{eqnarray}\begin{array}{lll}{\mu }_{\alpha }^{{\prime} } & = & -\displaystyle \frac{\partial U}{\partial {E}_{\alpha }}={\mu }_{\alpha }+{\alpha }_{\alpha \beta }^{{\prime} }{E}_{\beta }+\displaystyle \frac{1}{2}{a}_{\alpha ,\beta \gamma }^{{\prime} }{E}_{\beta \gamma }\\ & & +\,\displaystyle \frac{1}{6}{d}_{\alpha ,\beta \gamma \delta }^{{\prime} }{E}_{\beta \gamma \delta }+\cdots \end{array}\end{eqnarray}$(5b)$ \begin{eqnarray}\begin{array}{lll}{Q}_{\alpha \beta }^{{\prime} } & = & -2\displaystyle \frac{\partial U}{\partial {E}_{\alpha \beta }}={Q}_{\alpha \beta }+{a}_{\gamma ,\alpha \beta }^{{\prime} }{E}_{\gamma }\\ & & +\,\displaystyle \frac{1}{2}\left({c}_{\alpha \beta ,\gamma \delta }^{{\prime} }+{c}_{\gamma \delta ,\alpha \beta }^{{\prime} }\right){E}_{\gamma \delta }+\displaystyle \frac{1}{6}{g}_{\alpha \beta ,\gamma \delta \eta }^{{\prime} }{E}_{\gamma \delta \eta }+\,\cdots \\ & = & {Q}_{\alpha \beta }+{a}_{\gamma ,\alpha \beta }^{{\prime} }{E}_{\gamma }+{c}_{\alpha \beta ,\gamma \delta }^{{\prime} }{E}_{\gamma \delta }+\displaystyle \frac{1}{6}{g}_{\alpha \beta ,\gamma \delta \eta }^{{\prime} }{E}_{\gamma \delta \eta }+\,\cdots \end{array}\end{eqnarray}$(5c)$ \begin{eqnarray}\begin{array}{lll}{O}_{\alpha \beta \gamma }^{{\prime} } & = & -6\displaystyle \frac{\partial U}{\partial {E}_{\alpha \beta \gamma }}={O}_{\alpha \beta \gamma }+{d}_{\delta ,\alpha \beta \gamma }^{{\prime} }{E}_{\delta }+\displaystyle \frac{1}{2}{g}_{\delta \eta ,\alpha \beta \gamma }^{{\prime} }{E}_{\delta \eta }\\ & & +\,\displaystyle \frac{1}{6}\left({p}_{\alpha \beta \gamma ,\delta \eta \xi }^{{\prime} }+{p}_{\delta \eta \xi ,\alpha \beta \gamma }^{{\prime} }\right){E}_{\delta \eta \xi }+\,\cdots \\ & = & {O}_{\alpha \beta \gamma }+{d}_{\delta ,\alpha \beta \gamma }^{{\prime} }{E}_{\delta }+\displaystyle \frac{1}{2}{g}_{\delta \eta ,\alpha \beta \gamma }^{{\prime} }{E}_{\delta \eta }\\ & & +\,\displaystyle \frac{1}{3}{p}_{\alpha \beta \gamma ,\delta \eta \xi }^{{\prime} }{E}_{\delta \eta \xi }+\,\cdots \end{array}\end{eqnarray}$(5d)$ \begin{eqnarray}{H}_{\alpha \beta \gamma \delta }^{{\prime} }=-24\displaystyle \frac{\partial U}{\partial {E}_{\alpha \beta \gamma \delta }}={H}_{\alpha \beta \gamma \delta }+\,\cdots ,\end{eqnarray}$in which all but the first term on the RHS are the induced multipoles by the external electric field and gradients. α_αβ is dipole polarizability, a_γ,αβ is dipole–quadrupole polarizability, d_δ,αβγ is dipole–octopole polarizability, c_αβ,γδ =c_γδ,αβ is quadrupole–quadrupole polarizability, g_αβ,γδη is quadrupole–octopole polarizability and p_{αβγ,δηξ} =p_{δηξ,αβγ} is octopole–octopole polarizability, respectively, which are:(6a)$ \begin{eqnarray}{\alpha }_{\alpha \beta }=\displaystyle \frac{\partial {\mu }_{\alpha }}{\partial {E}_{\beta }}=-\displaystyle \frac{{\partial }^{2}U}{\partial {E}_{\alpha }\partial {E}_{\beta }}\end{eqnarray}$(6b)$ \begin{eqnarray}{a}_{\gamma ,\alpha \beta }=2\displaystyle \frac{\partial {\mu }_{\gamma }}{\partial {E}_{\alpha \beta }}=\displaystyle \frac{\partial {Q}_{\alpha \beta }}{\partial {E}_{\gamma }}=-2\displaystyle \frac{{\partial }^{2}U}{\partial {E}_{\gamma }\partial {E}_{\alpha \beta }}\end{eqnarray}$(6c)$ \begin{eqnarray}{d}_{\delta ,\alpha \beta \gamma }=6\displaystyle \frac{\partial {\mu }_{\delta }}{\partial {E}_{\alpha \beta \gamma }}=\displaystyle \frac{\partial {O}_{\alpha \beta \gamma }}{\partial {E}_{\delta }}=-6\displaystyle \frac{{\partial }^{2}U}{\partial {E}_{\delta }\partial {E}_{\alpha \beta \gamma }}\end{eqnarray}$(6d)$ \begin{eqnarray}{c}_{\alpha \beta ,\gamma \delta }=\displaystyle \frac{\partial {Q}_{\alpha \beta }}{\partial {E}_{\gamma \delta }}=-2\displaystyle \frac{{\partial }^{2}U}{\partial {E}_{\alpha \beta }\partial {E}_{\gamma \delta }}\end{eqnarray}$(6e)$ \begin{eqnarray}{g}_{\delta \eta ,\alpha \beta \gamma }=6\displaystyle \frac{\partial {Q}_{\delta \eta }}{\partial {E}_{\alpha \beta \gamma }}=2\displaystyle \frac{\partial {O}_{\alpha \beta \gamma }}{\partial {E}_{\delta \eta }}=-12\displaystyle \frac{{\partial }^{2}U}{\partial {E}_{\delta \eta }\partial {E}_{\alpha \beta \gamma }}\end{eqnarray}$(6f)$ \begin{eqnarray}{p}_{\alpha \beta \gamma ,\delta \eta \xi }=3\displaystyle \frac{\partial {O}_{\alpha \beta \gamma }}{\partial {E}_{\delta \eta \xi }}=-18\displaystyle \frac{{\partial }^{2}U}{\partial {E}_{\alpha \beta \gamma }\partial {E}_{\delta \eta \xi }},\end{eqnarray}$in which the multipole–multipole polarizabilities are all in the traced form. The primes on the superscripts of the (hyper)polarizabilities in equation (5) denote that they may all be inducible upon additional external perturbation. For example:(7a)$ \begin{eqnarray}{\alpha }_{\alpha \beta }^{{\prime} }={\alpha }_{\alpha \beta }+{\beta }_{\alpha \beta \gamma }^{{\prime} }{E}_{\gamma }+\displaystyle \frac{1}{2}{b}_{\alpha \beta ,\gamma \delta }^{{\prime} }{E}_{\gamma \delta }+\,\cdots \end{eqnarray}$(7b)$ \begin{eqnarray}{\beta }_{\alpha \beta \gamma }^{{\prime} }={\beta }_{\alpha \beta \gamma }+{\gamma }_{\alpha \beta \gamma \delta }^{{\prime} }{E}_{\delta }+\,\cdots \end{eqnarray}$(7c)$ \begin{eqnarray}\begin{array}{lll}{a}_{\gamma ,\alpha \beta }^{{\prime} } & = & {a}_{\gamma ,\alpha \beta }+{b}_{\gamma \delta ,\alpha \beta }^{{\prime} }{E}_{\delta }+\displaystyle \frac{1}{2}\left({c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)^{\prime} }+{c}_{\gamma ,\delta \eta ,\alpha \beta }^{(2)^{\prime} }\right){E}_{\delta \eta }+\,\cdots \\ & = & {a}_{\gamma ,\alpha \beta }+{b}_{\gamma \delta ,\alpha \beta }^{{\prime} }{E}_{\delta }+{c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)^{\prime} }{E}_{\delta \eta }+\,\cdots ,\end{array}\end{eqnarray}$in which β_αβγ and γ_αβγδ are the first and second dipole hyperpolarizability, respectively, b_αβ,γδ is dipole–dipole–quadrupole polarizability, and ${c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)}={c}_{\gamma ,\delta \eta ,\alpha \beta }^{(2)}$ is dipole–quadrupole–quadrupole polarizability, which are:(8a)$ \begin{eqnarray}{\beta }_{\alpha \beta \gamma }=\displaystyle \frac{\partial {\alpha }_{\alpha \beta }}{\partial {E}_{\gamma }}=-\displaystyle \frac{{\partial }^{3}U}{\partial {E}_{\alpha }\partial {E}_{\beta }\partial {E}_{\gamma }}\end{eqnarray}$(8b)$ \begin{eqnarray}{\gamma }_{\alpha \beta \gamma \delta }=\displaystyle \frac{\partial {\beta }_{\alpha \beta \gamma }}{\partial {E}_{\delta }}=-\displaystyle \frac{{\partial }^{4}U}{\partial {E}_{\alpha }\partial {E}_{\beta }\partial {E}_{\gamma }\partial {E}_{\delta }}\end{eqnarray}$(8c)$ \begin{eqnarray}{b}_{\gamma \delta ,\alpha \beta }=2\displaystyle \frac{\partial {\alpha }_{\gamma \delta }}{\partial {E}_{\alpha \beta }}=\displaystyle \frac{\partial {a}_{\gamma ,\alpha \beta }}{\partial {E}_{\delta }}=-2\displaystyle \frac{{\partial }^{3}U}{\partial {E}_{\gamma }\partial {E}_{\delta }\partial {E}_{\alpha \beta }}\end{eqnarray}$(8d)$ \begin{eqnarray}{c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)}=\displaystyle \frac{\partial {a}_{\gamma ,\alpha \beta }}{\partial {E}_{\delta \eta }}=-2\displaystyle \frac{{\partial }^{3}U}{\partial {E}_{\gamma }\partial {E}_{\alpha \beta }\partial {E}_{\delta \eta }}.\end{eqnarray}$In the Taylor series in equations (5) and (7), the Greek subscripts of the multipole–multipole polarizabilities are grouped, with different groups separated by a comma. The convention in this study is that the subscripts in the same group are permutable, while the exchange of the groups of the subscripts separated by a comma is not allowed. Thus, the singly grouped subscripts of electric fields, gradients, multipoles and dipole (hyper)polarizabilities are permutable. For the multiply grouped subscripts, taking dipole–quadrupole polarizability as example, a_γ,αβ and a_γ,βα are allowed in the summation of equations (5a)–(5b), but a_αβ,γ and a_βα,γ are not, although the four terms are all equal. Thus, such non-exchangeable group convention includes 27 a_γ,αβ terms in total, while allowing group exchange of a_αβ,γ merely doubles the number of terms, with the values of the terms of the former two times the corresponding terms of the latter. A similar argument also applies to b_γδ,αβ, d_δ,αβγ, g_δη,αβγ, etc. On the other hand, it is notable that the groups of the subscripts of c_αβ,γδ and p_{αβγ,δηξ}, with same number of subscripts in some groups, are inherently exchangeable. Thus, they are explicitly included in equations (5b)–(5c). The same argument also accounts for the ${c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)}$ in equation (7c), because the groups, αβ and δη, in the subscripts of ${c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)}$ are inherently exchangeable. Generally, if there are n inherently exchangeable groups in the subscripts of a multipole–multipole polarizability, the factor on the corresponding term is n! to account for the permutations of the inherently exchangeable groups.

It should be noted that the above convention was implicitly adopted by the Buckingham expansion in its originally proposed form [3–5], and we thus call it the Buckingham convention hereafter. Other conventions are also adopted in literatures [12–14]. Since the Buckingham expansion has been widely adopted in the community when reporting the multipole–multipole polarizabilities [8, 9, 11], other conventions are often converted to the Buckingham convention. One conversion between two different conventions is given in appendix A. Moreover, the Buckingham convention leads to a simple generalized Buckingham expansion in the traced form, as shown in appendix B.

Though the above Taylor series can be continued toward infinite terms, we hopefully include in equations (5)–(8) most of the important ones. The induced multipoles and (hyper)polarizabilities can then be obtained via successive integrations in a reversible manner [13]. For example, the induced component of the dipole polarizability α_αβ on external fields E_γ and E_δ, from equation (7b), is:(9)$ \begin{eqnarray}\begin{array}{lll}{\alpha }_{\alpha \beta }^{{\prime} }-{\alpha }_{\alpha \beta } & = & \displaystyle {\int }_{0}^{{E}_{\gamma }}{\beta }_{\alpha \beta \gamma }^{{\prime} }{\rm{d}}(\lambda {E}_{\gamma })\\ & = & \displaystyle {\int }_{0}^{1}\left({\beta }_{\alpha \beta \gamma }+{\gamma }_{\alpha \beta \gamma \delta }^{{\prime} }\lambda {E}_{\delta }+\cdots \right){E}_{\gamma }{\rm{d}}\lambda \\ & = & {\beta }_{\alpha \beta \gamma }{E}_{\gamma }+\displaystyle \frac{1}{2}{\gamma }_{\alpha \beta \gamma \delta }^{{\prime} }{E}_{\gamma }{E}_{\delta }+\,\cdots .\end{array}\end{eqnarray}$For the second dipole hyperpolarizability γ_αβγδ in the final expression in equation (9), the factor 1/2 arises from the integration over the parameter $\lambda \left(0\leqslant \lambda \leqslant 1\right)$ for a reversible polarization process [13]. Feeding the induction component of α_αβ in equation (9) back to equation (7a), we have:(7a′)$ \begin{eqnarray}{\alpha }_{\alpha \beta }^{{\prime} }={\alpha }_{\alpha \beta }+{\beta }_{\alpha \beta \gamma }{E}_{\gamma }+\displaystyle \frac{1}{2}{\gamma }_{\alpha \beta \gamma \delta }^{{\prime} }{E}_{\gamma }{E}_{\delta }+\displaystyle \frac{1}{2}{b}_{\alpha \beta ,\gamma \delta }^{{\prime} }{E}_{\gamma \delta }+\,\cdots ,\end{eqnarray}$which gives the differential dipole polarizability [6] α_αβ on external E_γ, E_δ, and E_γδ. Integrating equation (7a′) over E_β and equation (7c) over E_βγ with the reversible polarization process, respectively, and avoiding the double counting on the term involving dipole–dipole–quadrupole polarizability b_αβ,γδ, we can obtain the induced dipole μ^ind, i.e.:(10a)$ \begin{eqnarray}\begin{array}{lll}{\mu }_{\alpha }^{{\rm{ind}}} & = & \displaystyle {\int }_{0}^{1}\left({\alpha }_{\alpha \beta }+{\beta }_{\alpha \beta \gamma }\lambda {E}_{\gamma }+\displaystyle \frac{1}{2}{\gamma }_{\alpha \beta \gamma \delta }^{{\prime} }\lambda {E}_{\gamma }\lambda {E}_{\delta }+\,\cdots \right){E}_{\beta }{\rm{d}}\lambda \\ & & +\,\displaystyle \frac{1}{2}\displaystyle {\int }_{0}^{1}\left({a}_{\alpha ,\beta \gamma }+{b}_{\alpha \delta ,\beta \gamma }^{{\prime} }{E}_{\delta }+{c}_{\alpha ,\beta \gamma ,\delta \eta }^{(2)^{\prime} }\lambda {E}_{\delta \eta }+\,\cdots \right){E}_{\beta \gamma }{\rm{d}}\lambda \\ & = & {\alpha }_{\alpha \beta }{E}_{\beta }+\displaystyle \frac{1}{2}{\beta }_{\alpha \beta \gamma }{E}_{\beta }{E}_{\gamma }+\displaystyle \frac{1}{6}{\gamma }_{\alpha \beta \gamma \delta }^{{\prime} }{E}_{\beta }{E}_{\gamma }{E}_{\delta }+\,\cdots \\ & & +\,\displaystyle \frac{1}{2}{a}_{\alpha ,\beta \gamma }{E}_{\beta \gamma }+\displaystyle \frac{1}{2}{b}_{\alpha \delta ,\beta \gamma }^{{\prime} }{E}_{\delta }{E}_{\beta \gamma }\\ & & +\,\displaystyle \frac{1}{4}{c}_{\alpha ,\beta \gamma ,\delta \eta }^{(2)^{\prime} }{E}_{\beta \gamma }{E}_{\delta \eta }+\,\cdots \end{array}\end{eqnarray}$(10b)$ \begin{eqnarray}\begin{array}{lll}{Q}_{\alpha \beta }^{{\rm{ind}}} & = & \displaystyle {\int }_{0}^{1}\left({a}_{\gamma ,\alpha \beta }+{b}_{\gamma \delta ,\alpha \beta }^{{\prime} }\lambda {E}_{\delta }+{c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)^{\prime} }{E}_{\delta \eta }+\,\cdots \right){E}_{\gamma }{\rm{d}}\lambda \\ & = & {a}_{\gamma ,\alpha \beta }{E}_{\gamma }+\displaystyle \frac{1}{2}{b}_{\gamma \delta ,\alpha \beta }^{{\prime} }{E}_{\gamma }{E}_{\delta }+{c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)^{\prime} }{E}_{\gamma }{E}_{\delta \eta }+\,\cdots .\end{array}\end{eqnarray}$Feeding equation (10a) back to equations (5a), and (10b) to equation (5b), we have:(5a′)$ \begin{eqnarray}\begin{array}{lll}{\mu }_{\alpha }^{{\prime} } & = & {\mu }_{\alpha }+{\alpha }_{\alpha \beta }{E}_{\beta }+\displaystyle \frac{1}{2}{\beta }_{\alpha \beta \gamma }{E}_{\beta }{E}_{\gamma }\\ & & +\,\displaystyle \frac{1}{6}{\gamma }_{\alpha \beta \gamma \delta }^{{\prime} }{E}_{\beta }{E}_{\gamma }{E}_{\delta }+\,\cdots \\ & & +\,\displaystyle \frac{1}{2}{a}_{\alpha ,\beta \gamma }{E}_{\beta \gamma }+\displaystyle \frac{1}{2}{b}_{\alpha \delta ,\beta \gamma }^{{\prime} }{E}_{\delta }{E}_{\beta \gamma }\\ & & +\,\displaystyle \frac{1}{4}{c}_{\alpha ,\beta \gamma ,\delta \eta }^{(2)^{\prime} }{E}_{\beta \gamma }{E}_{\delta \eta }+\,\cdots \\ & & +\,\displaystyle \frac{1}{6}{d}_{\alpha ,\beta \gamma \delta }^{{\prime} }{E}_{\beta \gamma \delta }+\,\cdots \end{array}\end{eqnarray}$(5b′)$ \begin{eqnarray}\begin{array}{lll}{Q}_{\alpha \beta }^{{\prime} } & = & {Q}_{\alpha \beta }+{a}_{\gamma ,\alpha \beta }{E}_{\gamma }+\displaystyle \frac{1}{2}{b}_{\gamma \delta ,\alpha \beta }^{{\prime} }{E}_{\gamma }{E}_{\delta }\\ & & +\,{c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)^{\prime} }{E}_{\gamma }{E}_{\delta \eta }+\,\cdots \\ & & +\,{c}_{\alpha \beta ,\gamma \delta }^{{\prime} }{E}_{\gamma \delta }+\displaystyle \frac{1}{6}{g}_{\alpha \beta ,\gamma \delta \eta }^{{\prime} }{E}_{\gamma \delta \eta }+\,\cdots .\end{array}\end{eqnarray}$The Buckingham expansion in the traced form can be obtained via successive integration of equations (5a)–(5d), with equations (5a)–(5b) replaced by equations (5a′)–(5b′). Thus:(11a)$ \begin{eqnarray}\begin{array}{lll}U-{U}_{0} & = & -\displaystyle {\int }_{0}^{{\rm{1}}}{\mu }_{\alpha }^{{\prime} }{E}_{\alpha }{\rm{d}}\lambda =-\left({\mu }_{\alpha }+\displaystyle \frac{1}{2}{\alpha }_{\alpha \beta }{E}_{\beta }\right.\\ & & +\,\displaystyle \frac{1}{3!}{\beta }_{\alpha \beta \gamma }{E}_{\beta }{E}_{\gamma }+\displaystyle \frac{1}{4!}{\gamma }_{\alpha \beta \gamma \delta }{E}_{\beta }{E}_{\gamma }{E}_{\delta }\\ & & +\displaystyle \frac{1}{2}{a}_{\alpha ,\beta \gamma }{E}_{\beta \gamma }+\displaystyle \frac{1}{4}{b}_{\alpha \beta ,\gamma \delta }{E}_{\beta }{E}_{\gamma \delta }\\ & & \left.+\,\displaystyle \frac{1}{4}{c}_{\alpha ,\beta \gamma ,\delta \eta }^{(2)}{E}_{\beta \gamma }{E}_{\delta \eta }+\displaystyle \frac{1}{3!}{d}_{\alpha ,\beta \gamma \delta }{E}_{\beta \gamma \delta }+\cdots \right){E}_{\alpha }\end{array}\end{eqnarray}$(11b)$ \begin{eqnarray}\begin{array}{lll}U-{U}_{0} & = & -\displaystyle \frac{1}{2}\displaystyle {\int }_{0}^{1}{Q}_{\alpha \beta }^{{\prime} }{E}_{\alpha \beta }{\rm{d}}\lambda \\ & = & -\,\displaystyle \frac{1}{2}\left({Q}_{\alpha \beta }+{a}_{\gamma ,\alpha \beta }{E}_{\gamma }+\displaystyle \frac{1}{2}{b}_{\gamma \delta ,\alpha \beta }{E}_{\gamma }{E}_{\delta }\right.\\ & & +\,\displaystyle \frac{1}{2}{c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)}{E}_{\gamma }{E}_{\delta \eta }+\displaystyle \frac{1}{2}{c}_{\alpha \beta ,\gamma \delta }{E}_{\gamma \delta }\\ & & \left.+\,\displaystyle \frac{1}{3!}{g}_{\alpha \beta ,\gamma \delta \eta }{E}_{\gamma \delta \eta }+\,\cdots \right){E}_{\alpha \beta }\end{array}\end{eqnarray}$(11c)$ \begin{eqnarray}\begin{array}{lll}U-{U}_{0} & = & -\,\displaystyle \frac{1}{6}\displaystyle {\int }_{0}^{1}{O}_{\alpha \beta \gamma }^{{\prime} }{E}_{\alpha \beta \gamma }{\rm{d}}\lambda \\ & = & -\,\displaystyle \frac{1}{6}\left({O}_{\alpha \beta \gamma }+{d}_{\delta ,\alpha \beta \gamma }{E}_{\delta }+\displaystyle \frac{1}{2}{g}_{\delta \eta ,\alpha \beta \gamma }{E}_{\delta \eta }\right.\\ & & \left.+\,\displaystyle \frac{1}{6}{p}_{\alpha \beta \gamma ,\delta \eta \xi }{E}_{\delta \eta \xi }+\,\cdots \right){E}_{\alpha \beta \gamma }\end{array}\end{eqnarray}$(11d)$ \begin{eqnarray}\begin{array}{lll}U-{U}_{0} & = & -\displaystyle \frac{1}{24}\displaystyle {\int }_{0}^{1}{H}_{\alpha \beta \gamma \delta }^{{\prime} }{E}_{\alpha \beta \gamma \delta }{\rm{d}}\lambda \\ & = & -\displaystyle \frac{1}{24}\left({H}_{\alpha \beta \gamma \delta }+\cdots \right){E}_{\alpha \beta \gamma \delta }\end{array}\end{eqnarray}$
$ \begin{eqnarray*}\vdots \end{eqnarray*}$in which the primes on the superscripts are dropped in the final expression, because they are permanent ones without additional external electric fields or gradients other than those listed in the otherwise presented ones. With such understanding, the superscript primes are hereafter dropped. Using equations (11), equation (2) can be further expanded to include both multipoles and (hyper)polarizabilities. Summing equation (11) and skipping the same terms to avoid double count, we can expand the energy of a polarizable charge distribution. The resultant expansion is the Buckingham expansion in the traced form, i.e.:(12)$ \begin{eqnarray}\begin{array}{lll}U & = & {U}_{0}+q\varphi \\ & & -\,\left({\mu }_{\alpha }{E}_{\alpha }+\displaystyle \frac{1}{2}{\alpha }_{\alpha \beta }{E}_{\alpha }{E}_{\beta }+\displaystyle \frac{1}{6}{\beta }_{\alpha \beta \gamma }{E}_{\alpha }{E}_{\beta }{E}_{\gamma }\right.\\ & & \left.+\,\displaystyle \frac{1}{24}{\gamma }_{\alpha \beta \gamma \delta }{E}_{\alpha }{E}_{\beta }{E}_{\gamma }{E}_{\delta }+\,\cdots \right)\\ & & -\,\displaystyle \frac{1}{2}\left[\left({Q}_{\alpha \beta }+{a}_{\gamma ,\alpha \beta }{E}_{\gamma }+\displaystyle \frac{1}{2}{b}_{\gamma \delta ,\alpha \beta }{E}_{\gamma }{E}_{\delta }+\,\cdots \right)\right.\\ & & \left.+\,\displaystyle \frac{1}{2}\left({c}_{\alpha \beta ,\gamma \delta }+{c}_{\eta ,\alpha \beta ,\gamma \delta }^{(2)}{E}_{\eta }+\cdots \right){E}_{\gamma \delta }+\,\cdots \right]{E}_{\alpha \beta }\\ & & -\,\displaystyle \frac{1}{6}\left({O}_{\alpha \beta \gamma }+{d}_{\delta ,\alpha \beta \gamma }{E}_{\delta }+\displaystyle \frac{1}{2}{g}_{\delta \eta ,\alpha \beta \gamma }{E}_{\delta \eta }\right.\\ & & \left.+\,\displaystyle \frac{1}{6}{p}_{\delta \eta \xi ,\alpha \beta \gamma }{E}_{\delta \eta \xi }+\,\cdots \right){E}_{\alpha \beta \gamma }\\ & & -\,\displaystyle \frac{1}{24}{H}_{\alpha \beta \gamma \delta }{E}_{\alpha \beta \gamma \delta }-\,\cdots .\end{array}\end{eqnarray}$From equation (12), we can deduce that for a polarizable charge distribution the Buckingham expansion in the traced form can be constructed via successive Taylor series toward infinite terms, and all the multipoles in equation (4) and the multipole–multipole polarizabilities in equations (6) and (8) are the traced ones. Note that the Buckingham convention is adopted throughout the derivation. In appendix B, we give the general formula for the complete Buckingham expansion in the traced form. In the actual numerical calculations, it is convenient to get the traced multipoles and polarizabilities, and then transform them to their traceless counterparts if needed, as discussed below.

3. The Buckingham expansion in the traceless form

Without losing generality, we consider the externally perturbed φ is generated by a test positive unit electron charge e at R far from O on which the charge distribution ρ(r) is located. The electrostatic potential at r due to e at R is φ(r)=e/∣R − r∣, which, upon binomial expansion with respect to O, is:(13)$ \begin{eqnarray}\begin{array}{lll}\varphi ({\boldsymbol{r}}) & = & e\left(\displaystyle \frac{1}{R}+{r}_{\alpha }\displaystyle \frac{\partial {R}^{-1}}{\partial {R}_{\alpha }}-\displaystyle \frac{1}{2}{r}_{\alpha }{r}_{\beta }\displaystyle \frac{{\partial }^{2}{R}^{-1}}{\partial {R}_{\alpha }\partial {R}_{\beta }}\right.\\ & & +\,\displaystyle \frac{1}{3!}{r}_{\alpha }{r}_{\beta }{r}_{\gamma }\displaystyle \frac{{\partial }^{3}{R}^{-1}}{\partial {R}_{\alpha }\partial {R}_{\beta }\partial {R}_{\gamma }}\\ & & \left.-\,\displaystyle \frac{1}{4!}{r}_{\alpha }{r}_{\beta }{r}_{\gamma }{r}_{\delta }\displaystyle \frac{{\partial }^{4}{R}^{-1}}{\partial {R}_{\alpha }\partial {R}_{\beta }\partial {R}_{\gamma }\partial {R}_{\delta }}+\,\cdots \right).\end{array}\end{eqnarray}$Comparing equation (13) with equation (1), we have the external potential, electric field, and gradients at O, due to e at R, i.e.:(14a)$ \begin{eqnarray}\varphi (0)=\displaystyle \frac{e}{R}\end{eqnarray}$(14b)$ \begin{eqnarray}{E}_{\alpha }=-e\displaystyle \frac{\partial {R}^{-1}}{\partial {R}_{\alpha }}=e\displaystyle \frac{{R}_{\alpha }}{{R}^{3}}\end{eqnarray}$(14c)$ \begin{eqnarray}{E}_{\alpha \beta }=e\displaystyle \frac{{\partial }^{2}{R}^{-1}}{\partial {R}_{\alpha }\partial {R}_{\beta }}=e\displaystyle \frac{1}{{R}^{5}}\left(3{R}_{\alpha }{R}_{\beta }-{R}^{2}{\delta }_{\alpha \beta }\right)\end{eqnarray}$(14d)$ \begin{eqnarray}\begin{array}{lll}{E}_{\alpha \beta \gamma } & = & -e\displaystyle \frac{{\partial }^{3}{R}^{-1}}{\partial {R}_{\alpha }\partial {R}_{\beta }\partial {R}_{\gamma }}\\ & = & e\displaystyle \frac{1}{{R}^{7}}\left[15{R}_{\alpha }{R}_{\beta }{R}_{\gamma }-3{R}^{2}\right.\\ & & \left.\times \,\left({R}_{\alpha }{\delta }_{\beta \gamma }+{R}_{\beta }{\delta }_{\alpha \gamma }+{R}_{\gamma }{\delta }_{\alpha \beta }\right)\right]\end{array}\end{eqnarray}$(14e)$ \begin{eqnarray}\begin{array}{lll}{E}_{\alpha \beta \gamma \delta } & = & e\displaystyle \frac{{\partial }^{4}{R}^{-1}}{\partial {R}_{\alpha }\partial {R}_{\beta }\partial {R}_{\gamma }\partial {R}_{\delta }}\\ & = & e\displaystyle \frac{1}{{R}^{9}}\left[105{R}_{\alpha }{R}_{\beta }{R}_{\gamma }{R}_{\delta }\right.\\ & & -\,15{R}^{2}\left({R}_{\alpha }{R}_{\beta }{\delta }_{\gamma \delta }+{R}_{\alpha }{R}_{\gamma }{\delta }_{\beta \delta }+{R}_{\alpha }{R}_{\delta }{\delta }_{\beta \gamma }\right.\\ & & \left.+\,{R}_{\beta }{R}_{\gamma }{\delta }_{\alpha \delta }+{R}_{\beta }{R}_{\delta }{\delta }_{\alpha \gamma }+{R}_{\gamma }{R}_{\delta }{\delta }_{\alpha \beta }\right)\\ & & \left.+\,3{R}^{4}\left({\delta }_{\alpha \beta }{\delta }_{\gamma \delta }+{\delta }_{\alpha \gamma }{\delta }_{\beta \delta }+{\delta }_{\alpha \delta }{\delta }_{\beta \gamma }\right)\right]\end{array}\end{eqnarray}$
$ \begin{eqnarray*}\vdots \end{eqnarray*}$in which ${\delta }_{\alpha \beta }$ denotes the Kronecker symbol, which is 1 if α=β and is 0 otherwise. The subscripts α, β, γ, δ on the LHS of equation (14) refer to a specific Cartesian component of electric field or gradients.

From equation (14), it is evident that the electric field gradients are symmetric with respect to the permutations of the subscripts. Also, they are traceless, i.e. E_αα =E_αββ =E_αβγγ =0, etc, with repeated dummy variable summing over $\{x,y,z\},$ because O and R are spatially well separated and the Laplace equation is satisfied on any place other than R where e is located [2, 4, 12]. We next define auxiliary external electric field gradients in the traced form, F_αβ, F_αβγ, F_αβγδ, with the aid of equations (14c)–(14e), i.e.:(15a)$ \begin{eqnarray}{F}_{\alpha \beta }=e\displaystyle \frac{3{R}_{\alpha }{R}_{\beta }}{{R}^{5}}\end{eqnarray}$(15b)$ \begin{eqnarray}{F}_{\alpha \beta \gamma }=e\displaystyle \frac{15{R}_{\alpha }{R}_{\beta }{R}_{\gamma }}{{R}^{7}}\end{eqnarray}$(15c)$ \begin{eqnarray}{F}_{\alpha \beta \gamma \delta }=e\displaystyle \frac{105{R}_{\alpha }{R}_{\beta }{R}_{\gamma }{R}_{\delta }}{{R}^{9}}.\end{eqnarray}$Since electric field E_α is vector, it does not possess a ‘trace’. Using equation (15), the electric field gradients in equations (14c)–(14e) can be written as:(16a)$ \begin{eqnarray}{E}_{\alpha \beta }={F}_{\alpha \beta }-\displaystyle \frac{1}{3}{F}_{\mu \mu }{\delta }_{\alpha \beta }\end{eqnarray}$(16b)$ \begin{eqnarray}{E}_{\alpha \beta \gamma }={F}_{\alpha \beta \gamma }-\displaystyle \frac{1}{5}\left({F}_{\alpha \mu \mu }{\delta }_{\beta \gamma }+{F}_{\mu \beta \mu }{\delta }_{\alpha \gamma }+{F}_{\mu \mu \gamma }{\delta }_{\alpha \beta }\right)\end{eqnarray}$(16c)$ \begin{eqnarray}\begin{array}{lll}{E}_{\alpha \beta \gamma \delta } & = & {F}_{\alpha \beta \gamma \delta }-\displaystyle \frac{1}{7}\left({F}_{\alpha \beta \mu \mu }{\delta }_{\gamma \delta }+{F}_{\alpha \mu \gamma \mu }{\delta }_{\beta \delta }+{F}_{\alpha \mu \mu \delta }{\delta }_{\beta \gamma }\right.\\ & & \left.+{F}_{\mu \beta \gamma \mu }{\delta }_{\alpha \delta }+{F}_{\mu \beta \mu \delta }{\delta }_{\alpha \gamma }+{F}_{\mu \mu \gamma \delta }{\delta }_{\alpha \beta }\right)\\ & & +\,\displaystyle \frac{1}{35}\left({F}_{\mu \mu \nu \nu }{\delta }_{\alpha \beta }{\delta }_{\gamma \delta }+{F}_{\mu \nu \mu \nu }{\delta }_{\alpha \gamma }{\delta }_{\beta \delta }+{F}_{\mu \nu \nu \mu }{\delta }_{\alpha \delta }{\delta }_{\beta \gamma }\right),\end{array}\end{eqnarray}$in which the places of the subscripts α, β, γ, δ on the RHS are consistent with those on the LHS, which represents a specific Cartesian component of electric field or gradient. Equation (16) is derived by using the identity ${R}^{2}={R}_{\mu }{R}_{\mu }={R}_{x}^{2}+{R}_{y}^{2}+{R}_{z}^{2},$ since the subscript $\mu \in \{x,y,z\}$ is dummy, and similarly R⁴=R_μR_μR_νR_ν =R_μR_νR_μR_ν =R_μR_νR_νR_μ, in which the subscripts μ and ν are both dummy. It is easy to verify the traceless feature of electric field gradients by either using equations (14c)–(14e), or, equivalently, using equation (16), i.e.:(17a)$ \begin{eqnarray}{E}_{\alpha \alpha }=0\end{eqnarray}$(17b)$ \begin{eqnarray}{E}_{\alpha \alpha \gamma }={E}_{\alpha \beta \alpha }={E}_{\alpha \beta \beta }=0\end{eqnarray}$(17c)$ \begin{eqnarray}{E}_{\alpha \alpha \gamma \delta }={E}_{\alpha \beta \alpha \delta }={E}_{\alpha \beta \gamma \alpha }={E}_{\alpha \beta \beta \delta }={E}_{\alpha \beta \gamma \beta }={E}_{\alpha \beta \gamma \gamma }=0.\end{eqnarray}$For the above identities, equation (17a) is obtained via equation (16a), i.e. ${E}_{\alpha \alpha }={F}_{\alpha \alpha }-\tfrac{1}{3}{F}_{\eta \eta }{\delta }_{\alpha \alpha }=0,$ using the identity δ_αα for the repeated dummy subscript α. Using the first identity E_ααγ in equation (17b) as example, we have from equation (16b):
$ \begin{eqnarray*}\begin{array}{lll}{E}_{\alpha \alpha \gamma } & = & {F}_{\alpha \alpha \gamma }-\displaystyle \frac{1}{5}\left({F}_{\alpha \mu \mu }{\delta }_{\alpha \gamma }+{F}_{\mu \alpha \mu }{\delta }_{\alpha \gamma }+{F}_{\mu \mu \gamma }{\delta }_{\alpha \alpha }\right)\\ & = & {F}_{\alpha \alpha \gamma }-\displaystyle \frac{1}{5}\left({F}_{\gamma \mu \mu }+{F}_{\mu \gamma \mu }+3{F}_{\mu \mu \gamma }\right)=0,\end{array}\end{eqnarray*}$
in which we use the fact that the repeated subscript α is dummy, so that δ_αα =3 and α can be replaced by another dummy variable μ, i.e. F_ααγ =F_μμγ. Additionally, F_γμμ =F_μγμ =F_μμγ since F_αβγ is invariant upon the permutation of the subscripts. Similar arguments apply on the other identities in equations (17b)–(17c).

The traceless feature of the electric field gradients imposes the traceless condition on the corresponding multipoles. Taking, again, E_αβ as an example, and feeding E_αβ in equation (16a) to the corresponding traced quadrupole Q_αβ contribution to the energy in equation (2), we have $-\tfrac{1}{2}{Q}_{\alpha \beta }{E}_{\alpha \beta }$$=-\tfrac{1}{2}{Q}_{\alpha \beta }\left({F}_{\alpha \beta }-\tfrac{1}{3}{F}_{\mu \mu }{\delta }_{\alpha \beta }\right)$$=\,-\tfrac{1}{2}\left({Q}_{\alpha \beta }{F}_{\alpha \beta }-\tfrac{1}{3}{Q}_{\alpha \alpha }{F}_{\mu \mu }\right).$ Since all the subscripts are dummy, the last term can be written as ${Q}_{\alpha \alpha }{F}_{\mu \mu }={Q}_{\mu \mu }{F}_{\alpha \alpha }={Q}_{\mu \mu }{\delta }_{\alpha \beta }{F}_{\alpha \beta }.$ Thus:(18)$ \begin{eqnarray}\begin{array}{lll}-\displaystyle \frac{1}{2}{Q}_{\alpha \beta }{E}_{\alpha \beta } & = & -\displaystyle \frac{1}{2}\left({Q}_{\alpha \beta }{F}_{\alpha \beta }-\displaystyle \frac{1}{3}{Q}_{\mu \mu }{\delta }_{\alpha \beta }{F}_{\alpha \beta }\right)\\ & = & -\,\displaystyle \frac{1}{2}\left({Q}_{\alpha \beta }-\displaystyle \frac{1}{3}{Q}_{\mu \mu }{\delta }_{\alpha \beta }\right){F}_{\alpha \beta }=-\displaystyle \frac{1}{2}{\theta }_{\alpha \beta }{F}_{\alpha \beta }.\end{array}\end{eqnarray}$It is obvious that ${\theta }_{\alpha \beta }={Q}_{\alpha \beta }-\tfrac{1}{3}{Q}_{\mu \mu }{\delta }_{\alpha \beta }$ on the RHS of equation (18) is traceless, i.e. θ_αα =0, for the same reasoning in equation (17). Thus, θ_αβ can be defined as the traceless quadrupole. On the other hand, such definition is not unique, because the traceless feature is maintained upon multiplication of any factor f_Θ, i.e. f_Θθ_αβ. The Gaussian package [15] apparently adopts f_Θ=1, while in nuclear physics it adopts f_Θ=3 [16]. Here, we adopt Buckingham’s definition with f_Θ=3/2, so that the traceless quadrupole is defined to be [4]:(19)$ \begin{eqnarray}{{\rm{\Theta }}}_{\alpha \beta }={f}_{{\rm{\Theta }}}{\theta }_{\alpha \beta }=\displaystyle \frac{1}{2}\left(3{Q}_{\alpha \beta }-{Q}_{\mu \mu }{\delta }_{\alpha \beta }\right).\end{eqnarray}$Thus, the traced Q_αβ in equation (3c) can be transformed to the traceless Θ_αβ with the above equation. Using equation (19) with f_Θ=3/2, the energy $-\tfrac{1}{2}{Q}_{\alpha \beta }{E}_{\alpha \beta },$ in the presence of solely external electric field gradient, needs to be scaled by $\tfrac{1}{{f}_{{\rm{\Theta }}}}=\tfrac{2}{3}.$ Using equations (11b), (18) and (19), we have:(20)$ \begin{eqnarray}\begin{array}{lll}U-{U}_{0} & = & -\displaystyle \frac{1}{2}{Q}_{\alpha \beta }{E}_{\alpha \beta }=-\displaystyle \frac{1}{2}\displaystyle \frac{1}{{f}_{{\rm{\Theta }}}}{f}_{{\rm{\Theta }}}{\theta }_{\alpha \beta }{F}_{\alpha \beta }\\ & = & -\,\displaystyle \frac{1}{3}{{\rm{\Theta }}}_{\alpha \beta }{F}_{\alpha \beta }=-\displaystyle \frac{1}{3}{{\rm{\Theta }}}_{\alpha \beta }{E}_{\alpha \beta }.\end{array}\end{eqnarray}$The validation of last expression of equation (20) is due to the traceless feature of E_αβ and Θ_αβ. Using equation (16a) and identity Θ_αα =0, the energy in external electric field gradient is:
$ \begin{eqnarray*}\begin{array}{lll}-\displaystyle \frac{1}{3}{{\rm{\Theta }}}_{\alpha \beta }{E}_{\alpha \beta } & = & -\displaystyle \frac{1}{3}{{\rm{\Theta }}}_{\alpha \beta }\left({F}_{\alpha \beta }-\displaystyle \frac{1}{3}{F}_{\mu \mu }{\delta }_{\alpha \beta }\right)\\ & = & -\,\displaystyle \frac{1}{3}{{\rm{\Theta }}}_{\alpha \beta }{F}_{\alpha \beta }+\displaystyle \frac{1}{9}{{\rm{\Theta }}}_{\alpha \alpha }{F}_{\mu \mu }=-\displaystyle \frac{1}{3}{{\rm{\Theta }}}_{\alpha \beta }{F}_{\alpha \beta }.\end{array}\end{eqnarray*}$
The above relation validates the last expression of equation (20).

From equation (20), we can see that it is equivalent to use traced Q_αβ with traceless E_αβ in equation (14c) or equation (16a), or to use traceless Θ_αβ, with the auxiliary traced F_αβ in equation (15a). On the other hand, using both traceless multipole and traceless electric field gradient is redundant. With such understanding, from ${Q}_{\alpha \beta }=-2\tfrac{\partial U}{\partial {E}_{\alpha \beta }},$ we see that when the charge distribution ρ(r) senses a variation of the traceless E_αβ, it feeds back the traced Q_αβ in response. In contrast, equation (20) also tells that ${{\rm{\Theta }}}_{\alpha \beta }=-3\tfrac{\partial U}{\partial {F}_{\alpha \beta }}.$ Thus, when the charge distribution ρ(r) senses a variation of the traced F_αβ, it feeds back the traceless Θ_αβ in response. Using equation (16a), we have:(21a)$ \begin{eqnarray}{J}_{\gamma \delta ,\alpha \beta }=\displaystyle \frac{\partial {E}_{\gamma \delta }}{\partial {F}_{\alpha \beta }}=\displaystyle \frac{\partial {F}_{\gamma \delta }}{\partial {F}_{\alpha \beta }}-\displaystyle \frac{1}{3}\displaystyle \frac{\partial {F}_{\mu \mu }}{\partial {F}_{\alpha \beta }}{\delta }_{\gamma \delta }={\delta }_{\alpha \gamma }{\delta }_{\beta \delta }-\displaystyle \frac{1}{3}{\delta }_{\alpha \beta }{\delta }_{\gamma \delta },\end{eqnarray}$in which $\tfrac{\partial {F}_{\mu \mu }}{\partial {F}_{\alpha \beta }}={\delta }_{\mu \alpha }{\delta }_{\mu \beta }={\delta }_{\alpha \beta }$ is used. The traceless quadrupole, Θ_αβ, can then be calculated, equivalently, via:(22a)$ \begin{eqnarray}\begin{array}{lll}{{\rm{\Theta }}}_{\alpha \beta } & = & -{f}_{{\rm{\Theta }}}2\displaystyle \frac{\partial U}{\partial {F}_{\alpha \beta }}=-\displaystyle \frac{3}{2}2\displaystyle \frac{\partial U}{\partial {E}_{\gamma \delta }}\displaystyle \frac{\partial {E}_{\gamma \delta }}{\partial {F}_{\alpha \beta }}\\ & = & \displaystyle \frac{3}{2}{Q}_{\gamma \delta }\left({\delta }_{\alpha \gamma }{\delta }_{\beta \delta }-\displaystyle \frac{1}{3}{\delta }_{\alpha \beta }{\delta }_{\gamma \delta }\right)\\ & = & \displaystyle \frac{1}{2}\left(3{Q}_{\alpha \beta }-{Q}_{\mu \mu }{\delta }_{\alpha \beta }\right),\end{array}\end{eqnarray}$in which ${Q}_{\alpha \beta }=-2\tfrac{\partial U}{\partial {E}_{\alpha \beta }}$ as in equation (5b), and f_Θ=3/2 as suggested by Buckingham is adopted. Thus, equation (22a) gives the same traceless Θ_αβ in equation (19), as expected.

The differentiation of energy with respect to the auxiliary traced electric field gradient in equation (15a) makes it simple and consistent for the transformation from the traced multipoles and multipole–multipole polarizabilities to the corresponding traceless counterparts. From equations (16b)–(16c), we have:(21b)$ \begin{eqnarray}\begin{array}{lll}{J}_{\delta \eta \xi ,\alpha \beta \gamma } & = & \displaystyle \frac{\partial {E}_{\delta \eta \xi }}{\partial {F}_{\alpha \beta \gamma }}={\delta }_{\alpha \delta }{\delta }_{\beta \eta }{\delta }_{\gamma \xi }\\ & & -\,\displaystyle \frac{1}{5}\left({\delta }_{\delta \alpha }{\delta }_{\beta \gamma }{\delta }_{\eta \xi }+{\delta }_{\beta \eta }{\delta }_{\alpha \gamma }{\delta }_{\delta \xi }+{\delta }_{\gamma \xi }{\delta }_{\alpha \beta }{\delta }_{\delta \eta }\right)\end{array}\end{eqnarray}$(21c)$ \begin{eqnarray}\begin{array}{lll}{J}_{\eta \xi \mu \nu ,\alpha \beta \gamma \delta } & = & \displaystyle \frac{\partial {E}_{\eta \xi \mu \nu }}{\partial {F}_{\alpha \beta \gamma \delta }}={\delta }_{\alpha \eta }{\delta }_{\beta \xi }{\delta }_{\gamma \mu }{\delta }_{\delta \nu }\\ & & -\,\displaystyle \frac{1}{7}\left({\delta }_{\alpha \eta }{\delta }_{\beta \xi }{\delta }_{\gamma \delta }\right.{\delta }_{\mu \nu }+{\delta }_{\alpha \eta }{\delta }_{\beta \delta }{\delta }_{\gamma \mu }{\delta }_{\xi \nu }\\ & & +\,{\delta }_{\alpha \eta }{\delta }_{\beta \gamma }{\delta }_{\delta \nu }{\delta }_{\xi \mu }+{\delta }_{\alpha \delta }{\delta }_{\beta \xi }{\delta }_{\gamma \mu }{\delta }_{\eta \nu }\\ & & \left.+\,{\delta }_{\alpha \gamma }{\delta }_{\beta \xi }{\delta }_{\delta \nu }{\delta }_{\eta \mu }+{\delta }_{\alpha \beta }{\delta }_{\gamma \mu }{\delta }_{\delta \nu }{\delta }_{\eta \xi }\right)\\ & & +\,\displaystyle \frac{1}{35}\left({\delta }_{\alpha \beta }{\delta }_{\gamma \delta }{\delta }_{\eta \xi }{\delta }_{\mu \nu }+{\delta }_{\alpha \gamma }{\delta }_{\beta \delta }{\delta }_{\eta \mu }{\delta }_{\xi \nu }\right.\\ & & \left.+\,{\delta }_{\alpha \delta }{\delta }_{\beta \gamma }{\delta }_{\eta \nu }{\delta }_{\xi \mu }\right).\end{array}\end{eqnarray}$Similar to equation (22a) for transforming the traced Q_αβ to the traceless Θ_αβ, the traced octopole O_αβγ and hexadecapole H_αβγδ in equations (5c)–(5d) can be transformed to the traceless octopole Ω_αβγ and hexadecapole Φ_αβγδ with the aid of equations (21b)–(21c), i.e.:(22b)$ \begin{eqnarray}\begin{array}{lll}{{\rm{\Omega }}}_{\alpha \beta \gamma } & = & -{f}_{{\rm{\Omega }}}6\displaystyle \frac{\partial U}{\partial {F}_{\alpha \beta \gamma }}=-\displaystyle \frac{15}{6}6\displaystyle \frac{\partial U}{\partial {E}_{\delta \eta \xi }}\displaystyle \frac{\partial {E}_{\delta \eta \xi }}{\partial {F}_{\alpha \beta \gamma }}\\ & = & \displaystyle \frac{1}{6}\left[15{O}_{\alpha \beta \gamma }-3\left({O}_{\alpha \mu \mu }{\delta }_{\beta \gamma }+{O}_{\mu \beta \mu }{\delta }_{\alpha \gamma }\right.\right.\\ & & \left.\left.+\,{O}_{\mu \mu \gamma }{\delta }_{\alpha \beta }\right)\right]\end{array}\end{eqnarray}$(22c)$ \begin{eqnarray}\begin{array}{lll}{{\rm{\Phi }}}_{\alpha \beta \gamma \delta } & = & -{f}_{{\rm{\Phi }}}24\displaystyle \frac{\partial U}{\partial {F}_{\alpha \beta \gamma \delta }}=-\displaystyle \frac{105}{24}24\displaystyle \frac{\partial U}{\partial {E}_{\eta \xi \mu \nu }}\displaystyle \frac{\partial {E}_{\eta \xi \mu \nu }}{\partial {F}_{\alpha \beta \gamma \delta }}\\ & = & \displaystyle \frac{1}{24}\left[105{H}_{\alpha \beta \gamma \delta }-15\left({H}_{\alpha \beta \mu \mu }{\delta }_{\gamma \delta }+{H}_{\alpha \mu \gamma \mu }{\delta }_{\beta \delta }\right.\right.\\ & & \left.+\,{H}_{\alpha \mu \mu \delta }{\delta }_{\beta \gamma }+{H}_{\mu \beta \gamma \mu }{\delta }_{\alpha \delta }+{H}_{\mu \beta \mu \delta }{\delta }_{\alpha \gamma }+{H}_{\mu \mu \gamma \delta }{\delta }_{\alpha \beta }\right)\\ & & \left.+\,3\left({H}_{\mu \mu \nu \nu }{\delta }_{\alpha \beta }{\delta }_{\gamma \delta }+{H}_{\mu \nu \mu \nu }{\delta }_{\alpha \gamma }{\delta }_{\beta \delta }+{H}_{\mu \nu \nu \mu }{\delta }_{\alpha \delta }{\delta }_{\beta \gamma }\right)\right].\end{array}\end{eqnarray}$In equations (22a)–(22c), the Buckingham factor of the transformation from traced 2^M-moment to the traceless 2^M-moment is:(23)$ \begin{eqnarray}{f}_{M}=\displaystyle \frac{(2M-1)!!}{M!},\end{eqnarray}$in which ! denotes factorial and !! denotes double factorial with $(2M-1)!!=(2M-1)\cdot $$(2M-3)\cdots 3\cdot 1.$ Thus, ${f}_{{\rm{\Omega }}}=\tfrac{15}{6}$ and ${f}_{{\rm{\Phi }}}=\tfrac{105}{24}$ are used in equations (22b) and (22c), respectively. If ρ(r) is known, the traceless multipoles can be written explicitly using equations (3c)–(3e) and (22a)–(22c), i.e.:(24a)$ \begin{eqnarray}{{\rm{\Theta }}}_{\alpha \beta }=\displaystyle \frac{1}{2}\displaystyle {\int }_{V}{\rm{d}}v\rho ({\boldsymbol{r}})\left(3{r}_{\alpha }{r}_{\beta }-{r}^{2}{\delta }_{\alpha \beta }\right)\end{eqnarray}$(24b)$ \begin{eqnarray}\begin{array}{lll}{{\rm{\Omega }}}_{\alpha \beta \gamma } & = & \displaystyle \frac{1}{6}\displaystyle {\int }_{V}{\rm{d}}v\rho ({\boldsymbol{r}})\left[15{r}_{\alpha }{r}_{\beta }{r}_{\gamma }-3{r}^{2}\left({r}_{\alpha }{\delta }_{\beta \gamma }\right.\right.\\ & & \left.+\,\left.{r}_{\beta }{\delta }_{\alpha \gamma }+{r}_{\gamma }{\delta }_{\alpha \beta }\right)\right]\end{array}\end{eqnarray}$(24c)$ \begin{eqnarray}\begin{array}{lll}{{\rm{\Phi }}}_{\alpha \beta \gamma \delta } & = & \displaystyle \frac{1}{24}\displaystyle {\int }_{V}{\rm{d}}v\rho ({\boldsymbol{r}})\left[105{r}_{\alpha }{r}_{\beta }{r}_{\gamma }{r}_{\delta }-15{r}^{2}\left({r}_{\alpha }{r}_{\beta }{\delta }_{\gamma \delta }\right.\right.\\ & & +\,{r}_{\alpha }{r}_{\gamma }{\delta }_{\beta \delta }+{r}_{\alpha }{r}_{\delta }{\delta }_{\beta \gamma }+{r}_{\beta }{r}_{\gamma }{\delta }_{\alpha \delta }+{r}_{\beta }{r}_{\delta }{\delta }_{\alpha \gamma }\\ & & \left.\left.+\,{r}_{\gamma }{r}_{\delta }{\delta }_{\alpha \beta }\right)+3{r}^{4}\left({\delta }_{\alpha \beta }{\delta }_{\gamma \delta }+{\delta }_{\alpha \gamma }{\delta }_{\beta \delta }+{\delta }_{\alpha \delta }{\delta }_{\beta \gamma }\right)\right].\end{array}\end{eqnarray}$Using the same reasoning as in equation (17), it can be verified by using either equation (22) or equation (24) for the following traceless identities:(25a)$ \begin{eqnarray}{{\rm{\Theta }}}_{\alpha \alpha }=0\end{eqnarray}$(25b)$ \begin{eqnarray}{{\rm{\Omega }}}_{\alpha \alpha \gamma }={{\rm{\Omega }}}_{\alpha \beta \alpha }={{\rm{\Omega }}}_{\alpha \beta \beta }=0\end{eqnarray}$(25c)$ \begin{eqnarray}{{\rm{\Phi }}}_{\alpha \alpha \gamma \delta }={{\rm{\Phi }}}_{\alpha \beta \alpha \delta }={{\rm{\Phi }}}_{\alpha \beta \gamma \alpha }={{\rm{\Phi }}}_{\alpha \beta \beta \delta }={{\rm{\Phi }}}_{\alpha \beta \gamma \beta }={{\rm{\Phi }}}_{\alpha \beta \gamma \gamma }=0.\end{eqnarray}$Though the transformation from the traced to the traceless multipoles is somehow tedious, the expression of energy is straightforward, as shown in equation (20) for the quadrupole. It can be seen that the energy would be the same by multiplying and then dividing f_M in equation (23). Extending the equation (20), which is for the permanent quadrupole, to the corresponding inducible counterpart in equation (11b), we have:(26a)$ \begin{eqnarray}\begin{array}{lll}U-{U}^{0} & = & -\displaystyle \frac{1}{{f}_{{\rm{\Theta }}}}\displaystyle \frac{1}{2}\left({f}_{{\rm{\Theta }}}{Q}_{\alpha \beta }+{f}_{{\rm{\Theta }}}{a}_{\gamma ,\alpha \beta }{E}_{\gamma }+\displaystyle \frac{1}{2}{f}_{{\rm{\Theta }}}{b}_{\gamma \delta ,\alpha \beta }{E}_{\gamma }{E}_{\delta }\right.\\ & & +\,\displaystyle \frac{1}{2}{f}_{{\rm{\Theta }}}{c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)}{E}_{\gamma }{E}_{\delta \eta }+\displaystyle \frac{1}{2}{f}_{{\rm{\Theta }}}{c}_{\alpha \beta ,\gamma \delta }{E}_{\gamma \delta }\\ & & \left.+\,\displaystyle \frac{1}{{f}_{{\rm{\Omega }}}}\displaystyle \frac{1}{3!}{f}_{{\rm{\Theta }}}{f}_{{\rm{\Omega }}}{g}_{\alpha \beta ,\gamma \delta \eta }{E}_{\gamma \delta \eta }+\cdots \right){E}_{\alpha \beta }\\ & = & -\,\displaystyle \frac{1}{3}\left({{\rm{\Theta }}}_{\alpha \beta }+{A}_{\gamma ,\alpha \beta }{E}_{\gamma }+\displaystyle \frac{1}{2}{B}_{\gamma \delta ,\alpha \beta }{E}_{\gamma }{E}_{\delta }\right.\\ & & +\,\displaystyle \frac{1}{2}{C}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)}{E}_{\gamma }{E}_{\delta \eta }+\displaystyle \frac{1}{2}{C}_{\alpha \beta ,\gamma \delta }{E}_{\gamma \delta }\\ & & \left.+\,\displaystyle \frac{1}{15}{G}_{\alpha \beta ,\gamma \delta \eta }{E}_{\gamma \delta \eta }+\cdots \right){E}_{\alpha \beta },\end{array}\end{eqnarray}$in which ${f}_{{\rm{\Theta }}}=\tfrac{3}{2}$ and ${f}_{{\rm{\Omega }}}=\tfrac{15}{6},$ according to equation (23), and A_γ,αβ, B_γδ,αβ, C_αβ,γδ, ${C}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)},$G_αβ,γδη denote the traceless multipole–multipole polarizabilities, with the corresponding lowercase symbols denote the traced counterparts. Thus, the way to convert the energy from the expression of the traced multipoles and (hyper)polarizabilities to the expression of the corresponding traceless counterpart is straightforward, i.e. by multiplying f_Θ if electric field gradient E_αβ is involved, and multiplying f_Ω if electric field gradient E_αβγ is involved, etc. Such manipulation converts equation (11b) for the traced multipoles and multipole–multipole polarizabilities to the traceless counterparts in equation (26a). It is notable that for the terms involving quadrupole–quadrupole polarizability c_αβ,γδ and dipole–quadrupole–quadrupole polarizability ${c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)},$ which have an inherently exchangeable group in the sub-indices, the corresponding f_Θ is multiplied only once. This is because the corresponding traceless ${{\rm{\Theta }}}_{\alpha \beta }\,{\rm{a}}{\rm{n}}{\rm{d}}\,{A}_{\gamma ,\alpha \beta }$ energy expressions are obtained by a single multiplication with ${E}_{\gamma \delta }$_αβ, and both are converted to the traceless form, as discussed in equations (18)–(20). Moreover, it should be emphasized that equation (26a) merely denotes $-\tfrac{1}{2}{Q}_{\alpha \beta }^{{\prime} }{E}_{\alpha \beta }=-\tfrac{1}{3}{{\rm{\Theta }}}_{\alpha \beta }^{{\prime} }{E}_{\alpha \beta }$ by extending equation (20) to include inducible quadrupoles contributed from various external electric fields and gradients, but cannot give the explicit form of the multipole–multipole polarizabilities. For example, equation (26a) denotes $-\tfrac{1}{4}{c}_{\gamma \delta ,\alpha \beta }{E}_{\alpha \beta }\,=-\tfrac{1}{6}{C}_{\gamma \delta ,\alpha \beta }{E}_{\alpha \beta },$ but cannot give the explicit form of C_γδ,αβ, which will be derived below in equation (28).

Based on the above discussions regarding equation (26a), the contribution of the energy from the traceless octopole and hexadecapole, with inducible components can be extended from equations (11c)–(11d), i.e.:(26b)$ \begin{eqnarray}\begin{array}{lll}U-{U}^{0} & = & -\displaystyle \frac{1}{{f}_{{\rm{\Omega }}}}\displaystyle \frac{1}{6}\left({f}_{{\rm{\Omega }}}{O}_{\alpha \beta \gamma }+{f}_{{\rm{\Omega }}}{d}_{\delta ,\alpha \beta \gamma }{E}_{\delta }\right.\\ & & +\,\displaystyle \frac{1}{{f}_{{\rm{\Theta }}}}\displaystyle \frac{1}{2}{f}_{{\rm{\Theta }}}{f}_{{\rm{\Omega }}}{g}_{\delta \eta ,\alpha \beta \gamma }{E}_{\delta \eta }\\ & & \left.+\,\displaystyle \frac{1}{3!}{f}_{{\rm{\Omega }}}{p}_{\alpha \beta \gamma ,\delta \eta \xi }{E}_{\delta \eta \xi }+\cdots \right){E}_{\alpha \beta \gamma }\\ & = & -\,\displaystyle \frac{1}{15}\left({{\rm{\Omega }}}_{\alpha \beta \gamma }+{D}_{\delta ,\alpha \beta \gamma }{E}_{\delta }\right.\\ & & \left.+\,\displaystyle \frac{1}{{\rm{3}}}{G}_{\delta \eta ,\alpha \beta \gamma }{E}_{\delta \eta }+\displaystyle \frac{1}{6}{P}_{\alpha \beta \gamma ,\delta \eta \xi }{E}_{\delta \eta \xi }+\cdots \right)\\ & & \times \,{E}_{\alpha \beta \gamma }\end{array}\end{eqnarray}$(26c)$ \begin{eqnarray}\begin{array}{lll}U-{U}^{0} & = & -\displaystyle \frac{1}{{f}_{{\rm{\Phi }}}}\displaystyle \frac{1}{24}\left({f}_{{\rm{\Phi }}}{H}_{\alpha \beta \gamma \delta }+\cdots \right){E}_{\alpha \beta \gamma \delta }\\ & = & -\,\displaystyle \frac{1}{105}\left({{\rm{\Phi }}}_{\alpha \beta \gamma \delta }+\cdots \right){E}_{\alpha \beta \gamma \delta },\end{array}\end{eqnarray}$in which D_δ,αβγ and P_{δηξ,αβγ} are the dipole–octopole and octopole–octopole polarizabilities, respectively, in the traceless form. Note that since the groups of the subscripts of octopole–octopole polarizabilities are inherently exchangeable, the term involving it is multiplied by f_Ω only once. Summing equations (26a)–(26c), and skipping the same terms to avoid double counts, the Buckingham expansion of the electrostatic energy of a polarizable charge distribution in the traced form, equation (12), can then be transform to the traceless form, i.e.:(27)$ \begin{eqnarray}\begin{array}{lll}U & = & {U}^{0}+q\varphi \\ & & -\,\left({\mu }_{\alpha }{E}_{\alpha }+\displaystyle \frac{1}{2}{\alpha }_{\alpha \beta }{E}_{\alpha }{E}_{\beta }+\displaystyle \frac{1}{6}{\beta }_{\alpha \beta \gamma }{E}_{\alpha }{E}_{\beta }{E}_{\gamma }\right.\\ & & \left.+\,\displaystyle \frac{1}{24}{\gamma }_{\alpha \beta \gamma \delta }{E}_{\alpha }{E}_{\beta }{E}_{\gamma }{E}_{\delta }+\,\cdots \right)\\ & & -\,\displaystyle \frac{1}{3}\left[\left({{\rm{\Theta }}}_{\alpha \beta }+{A}_{\gamma ,\alpha \beta }{E}_{\gamma }+\displaystyle \frac{1}{2}{B}_{\gamma \delta ,\alpha \beta }{E}_{\gamma }{E}_{\delta }+\,\cdots \right)\right.\\ & & \left.+\,\displaystyle \frac{1}{2}\left({C}_{\alpha \beta ,\delta \eta }+{C}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)}{E}_{\gamma }+\,\cdots \right){E}_{\delta \eta }+\,\cdots \right]{E}_{\alpha \beta }\\ & & -\,\displaystyle \frac{1}{15}\left({{\rm{\Omega }}}_{\alpha \beta \gamma }+{D}_{\delta ,\alpha \beta \gamma }{E}_{\delta }+\displaystyle \frac{1}{{\rm{3}}}{G}_{\delta \eta ,\alpha \beta \gamma }{E}_{\delta \eta }\right.\\ & & \left.+\,\displaystyle \frac{1}{6}{P}_{\delta \eta \xi ,\alpha \beta \gamma }{E}_{\delta \eta \xi }+\,\cdots \right){E}_{\alpha \beta \gamma }\\ & & -\,\displaystyle \frac{1}{105}{{\rm{\Phi }}}_{\alpha \beta \gamma \delta }{E}_{\alpha \beta \gamma \delta }-\,\cdots .\end{array}\end{eqnarray}$Equation (27) is the final expression of the Buckingham expansion in the traceless form [3–5], in which all the multipoles and multipole–multipole polarizabilities are traceless.

The transformation from the traced multipoles to the traceless multipoles are given by equations (22a)–(22c). We next derive the transformation from the traced multipole–multipole polarizabilities in equations (6) and (8) to the traceless counterparts, using the traceless multipoles and the auxiliary traced electric field gradients. Using equation (27), we have:(28a)$ \begin{eqnarray}\begin{array}{lll}{A}_{\gamma ,\alpha \beta } & = & -3\displaystyle \frac{{\partial }^{2}U}{\partial {F}_{\alpha \beta }\partial {E}_{\gamma }}=\displaystyle \frac{\partial {{\rm{\Theta }}}_{\alpha \beta }}{\partial {E}_{\gamma }}\\ & = & \displaystyle \frac{1}{2}\left(3\displaystyle \frac{\partial {Q}_{\alpha \beta }}{\partial {E}_{\gamma }}-\displaystyle \frac{\partial {Q}_{\mu \mu }}{\partial {E}_{\gamma }}{\delta }_{\alpha \beta }\right)\\ & = & \displaystyle \frac{1}{2}\left(3{a}_{\gamma ,\alpha \beta }-{a}_{\gamma ,\mu \mu }{\delta }_{\alpha \beta }\right)\end{array}\end{eqnarray}$(28b)$ \begin{eqnarray}\begin{array}{lll}{B}_{\gamma \delta ,\alpha \beta } & = & -3\displaystyle \frac{{\partial }^{3}U}{\partial {F}_{\alpha \beta }\partial {E}_{\gamma }\partial {E}_{\delta }}\\ & = & \displaystyle \frac{\partial {A}_{\gamma ,\alpha \beta }}{\partial {E}_{\delta }}\\ & = & \displaystyle \frac{1}{2}\left(3{b}_{\gamma \delta ,\alpha \beta }-{b}_{\gamma \delta ,\mu \mu }{\delta }_{\alpha \beta }\right)\end{array}\end{eqnarray}$(28c)$ \begin{eqnarray}\begin{array}{lll}{C}_{\alpha \beta ,\gamma \delta } & = & -3\displaystyle \frac{{\partial }^{2}U}{\partial {F}_{\alpha \beta }\partial {F}_{\gamma \delta }}=\displaystyle \frac{\partial {{\rm{\Theta }}}_{\alpha \beta }}{\partial {F}_{\gamma \delta }}=\displaystyle \frac{\partial {{\rm{\Theta }}}_{\alpha \beta }}{\partial {E}_{\kappa \lambda }}\displaystyle \frac{\partial {E}_{\kappa \lambda }}{\partial {F}_{\gamma \delta }}\\ & = & \displaystyle \frac{1}{2}\left(3\displaystyle \frac{\partial {Q}_{\alpha \beta }}{\partial {E}_{\kappa \lambda }}-\displaystyle \frac{\partial {Q}_{\mu \mu }}{\partial {E}_{\kappa \lambda }}{\delta }_{\alpha \beta }\right){J}_{\kappa \lambda ,\gamma \delta }\\ & = & \displaystyle \frac{1}{6}\left(9{c}_{\alpha \beta ,\gamma \delta }-3{c}_{\alpha \beta ,\mu \mu }{\delta }_{\gamma \delta }-3{c}_{\mu \mu ,\gamma \delta }{\delta }_{\alpha \beta }\right.\\ & & \left.+\,{c}_{\mu \mu ,\nu \nu }{\delta }_{\alpha \beta }{\delta }_{\gamma \delta }\right)\end{array}\end{eqnarray}$(28d)$ \begin{eqnarray}\begin{array}{lll}{C}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)} & = & -3\displaystyle \frac{{\partial }^{3}U}{\partial {F}_{\alpha \beta }\partial {F}_{\delta \eta }\partial {E}_{\gamma }}=\displaystyle \frac{\partial {C}_{\alpha \beta ,\delta \eta }}{\partial {E}_{\gamma }}\\ & = & \displaystyle \frac{1}{6}\left(9{c}_{\gamma ,\alpha \beta ,\delta \eta }^{(2)}-3{c}_{\gamma ,\alpha \beta ,\mu \mu }^{(2)}{\delta }_{\delta \eta }\right.\\ & & \left.-\,3{c}_{\gamma ,\mu \mu ,\delta \eta }^{(2)}{\delta }_{\alpha \beta }+{c}_{\gamma ,\mu \mu ,\nu \nu }^{(2)}{\delta }_{\alpha \beta }{\delta }_{\delta \eta }\right)\end{array}\end{eqnarray}$(28e)$ \begin{eqnarray}\begin{array}{lll}{D}_{\delta ,\alpha \beta \gamma } & = & -15\displaystyle \frac{{\partial }^{2}U}{\partial {F}_{\alpha \beta \gamma }\partial {E}_{\delta }}=\displaystyle \frac{\partial {{\rm{\Omega }}}_{\alpha \beta \gamma }}{\partial {E}_{\delta }}\\ & = & \displaystyle \frac{1}{6}\left[15{d}_{\delta ,\alpha \beta \gamma }-3\left({d}_{\delta ,\alpha \mu \mu }{\delta }_{\beta \gamma }+{d}_{\delta ,\mu \beta \mu }{\delta }_{\alpha \gamma }\right.\right.\\ & & \left.+\,\left.{d}_{\delta ,\mu \mu \gamma }{\delta }_{\alpha \beta }\right)\right]\end{array}\end{eqnarray}$(28f)$ \begin{eqnarray}\begin{array}{lll}{G}_{\delta \eta ,\alpha \beta \gamma } & = & -{\rm{45}}\displaystyle \frac{{\partial }^{2}U}{\partial {F}_{\alpha \beta \gamma }\partial {F}_{\delta \eta }}={\rm{3}}\displaystyle \frac{\partial {{\rm{\Omega }}}_{\alpha \beta \gamma }}{\partial {F}_{\delta \eta }}={\rm{3}}\displaystyle \frac{\partial {{\rm{\Omega }}}_{\alpha \beta \gamma }}{\partial {E}_{\kappa \lambda }}{J}_{\kappa \lambda ,\delta \eta }\\ & = & \displaystyle \frac{1}{{\rm{4}}}\left[15{g}_{\kappa \lambda ,\alpha \beta \gamma }-3\left({g}_{\kappa \lambda ,\alpha \mu \mu }{\delta }_{\beta \gamma }+{g}_{\kappa \lambda ,\mu \beta \mu }{\delta }_{\alpha \gamma }\right.\right.\\ & & \left.\left.+\,{g}_{\kappa \lambda ,\mu \mu \gamma }{\delta }_{\alpha \beta }\right)\right]\left({\delta }_{\delta \kappa }{\delta }_{\eta \lambda }-\displaystyle \frac{1}{3}{\delta }_{\delta \eta }{\delta }_{\kappa \lambda }\right)\\ & = & \displaystyle \frac{1}{{\rm{4}}}\left[15{g}_{\delta \eta ,\alpha \beta \gamma }-5{g}_{\mu \mu ,\alpha \beta \gamma }{\delta }_{\delta \eta }-3\right.\\ & & \times \,\left({g}_{\delta \eta ,\alpha \nu \nu }{\delta }_{\beta \gamma }+{g}_{\delta \eta ,\nu \beta \nu }{\delta }_{\alpha \gamma }+{g}_{\delta \eta ,\nu \nu \gamma }{\delta }_{\alpha \beta }\right)\\ & & +\,{g}_{\mu \mu ,\alpha \nu \nu }{\delta }_{\beta \gamma }{\delta }_{\delta \eta }+{g}_{\mu \mu ,\nu \beta \nu }{\delta }_{\alpha \gamma }{\delta }_{\delta \eta }\\ & & \left.+\,{g}_{\mu \mu ,\nu \nu \gamma }{\delta }_{\alpha \beta }{\delta }_{\delta \eta }\right]\end{array}\end{eqnarray}$(28g)$ \begin{eqnarray}\begin{array}{lll}{P}_{\alpha \beta \gamma ,\delta \eta \xi } & = & -45\displaystyle \frac{{{\rm{\partial }}}^{2}U}{{\rm{\partial }}{F}_{\alpha \beta \gamma }{\rm{\partial }}{F}_{\delta \eta \xi }}=3\displaystyle \frac{{\rm{\partial }}{{\rm{\Omega }}}_{\alpha \beta \gamma }}{{\rm{\partial }}{F}_{\delta \eta \xi }}=3\displaystyle \frac{{\rm{\partial }}{{\rm{\Omega }}}_{\alpha \beta \gamma }}{{\rm{\partial }}{E}_{\lambda \mu \nu }}{J}_{\lambda \mu \nu ,\delta \eta \xi }\\ & = & \displaystyle \frac{1}{6}\left[15{p}_{\alpha \beta \gamma ,\delta \eta \xi }-3\left({p}_{\alpha \beta \gamma ,\delta \mu \mu }{\delta }_{\eta \xi }+{p}_{\alpha \beta \gamma ,\mu \eta \mu }{\delta }_{\delta \xi }\right.\right.\\ & & +\,{p}_{\alpha \beta \gamma ,\mu \mu \xi }{\delta }_{\delta \eta }+{p}_{\alpha \mu \mu ,\delta \eta \xi }{\delta }_{\beta \gamma }+{p}_{\mu \beta \mu ,\delta \eta \xi }{\delta }_{\alpha \gamma }\\ & & \left.+\,{p}_{\mu \mu \gamma ,\delta \eta \xi }{\delta }_{\alpha \beta }\right)+\displaystyle \frac{3}{5}\left({p}_{\alpha \mu \mu ,\delta \nu \nu }{\delta }_{\beta \gamma }{\delta }_{\eta \xi }+{p}_{\mu \beta \mu ,\delta \nu \nu }{\delta }_{\alpha \gamma }{\delta }_{\eta \xi }\right.\\ & & +\,{p}_{\mu \mu \gamma ,\delta \nu \nu }{\delta }_{\alpha \beta }{\delta }_{\eta \xi }+{p}_{\alpha \mu \mu ,\nu \eta \nu }{\delta }_{\beta \gamma }{\delta }_{\delta \xi }\\ & & +\,{p}_{\mu \beta \mu ,\nu \eta \nu }{\delta }_{\alpha \gamma }{\delta }_{\delta \xi }+{p}_{\mu \mu \gamma ,\nu \eta \nu }{\delta }_{\alpha \beta }{\delta }_{\delta \xi }+{p}_{\alpha \mu \mu ,\nu \nu \xi }{\delta }_{\beta \gamma }{\delta }_{\delta \eta }\\ & & \left.\left.+\,{p}_{\mu \beta \mu ,\nu \nu \xi }{\delta }_{\alpha \gamma }{\delta }_{\delta \eta }+{p}_{\mu \mu \gamma ,\nu \nu \xi }{\delta }_{\alpha \beta }{\delta }_{\delta \eta }\right)\right],\end{array}\end{eqnarray}$in which the chain rule is used in order to get the traceless multipole–multipole polarizabilities, with the aid of equations (21a)–(21b). Note that some of the transformations in equation (28), for example, equation (28c), have also been given by Applequist [13] using a different approach. It can be verified using equation (28) by the following traceless identities:(29a)$ \begin{eqnarray}{A}_{\gamma ,\alpha \alpha }=0\end{eqnarray}$(29b)$ \begin{eqnarray}{B}_{\gamma \delta ,\alpha \alpha }=0\end{eqnarray}$(29c)$ \begin{eqnarray}{C}_{\alpha \alpha ,\gamma \delta }={C}_{\alpha \beta ,\gamma \gamma }=0\end{eqnarray}$(29d)$ \begin{eqnarray}{C}_{\gamma ,\alpha \alpha ,\delta \eta }^{(2)}={C}_{\gamma ,\alpha \beta ,\delta \delta }^{(2)}=0\end{eqnarray}$(29e)$ \begin{eqnarray}{D}_{\delta ,\alpha \alpha \gamma }={D}_{\delta ,\alpha \beta \alpha }={D}_{\delta ,\alpha \beta \beta }=0\end{eqnarray}$(29f)$ \begin{eqnarray}{G}_{\delta \eta ,\alpha \alpha \gamma }={G}_{\delta \eta ,\alpha \beta \alpha }={G}_{\delta \eta ,\alpha \beta \beta }={G}_{\delta \delta ,\alpha \beta \gamma }=0\end{eqnarray}$(29g)$ \begin{eqnarray}\begin{array}{lll}{P}_{\alpha \alpha \gamma ,\delta \eta \xi } & = & {P}_{\alpha \beta \alpha ,\delta \eta \xi }={P}_{\alpha \beta \beta ,\delta \eta \xi }={P}_{\alpha \beta \gamma ,\delta \delta \xi }\\ & = & {P}_{\alpha \beta \gamma ,\delta \eta \delta }={P}_{\alpha \beta \gamma ,\delta \eta \eta }=0.\end{array}\end{eqnarray}$In obtaining the above identities, using the first identity D_δ,ααγ in equation (29e) as an example, we have from equation (28e):
$ \begin{eqnarray*}\begin{array}{lll}{D}_{\delta ,\alpha \alpha \gamma } & = & \displaystyle \frac{1}{6}\left[15{d}_{\delta ,\alpha \alpha \gamma }-3\left({d}_{\delta ,\alpha \mu \mu }{\delta }_{\alpha \gamma }+{d}_{\delta ,\mu \alpha \mu }{\delta }_{\alpha \gamma }\right.\right.\\ & & \left.+\,\left.{d}_{\delta ,\mu \mu \gamma }{\delta }_{\alpha \alpha }\right)\right]\\ & = & \displaystyle \frac{1}{6}\left[15{d}_{\delta ,\mu \mu \gamma }-3\left({d}_{\delta ,\gamma \mu \mu }+{d}_{\delta ,\mu \gamma \mu }+3{d}_{\delta ,\mu \mu \gamma }\right)\right]\\ & = & 0,\end{array}\end{eqnarray*}$
in which we use the fact that the repeated sub-index α is dummy, so that δ_αα =3 and α can be replaced by another dummy variable μ, i.e. D_δ,ααγ =D_δ,μμγ; and D_δ,γμμ =D_δ,μγμ =D_δ,μμγ since the multipole–multipole polarizabilities are invariant upon the permutation of the sub-indices in the same group. Similar arguments apply on the identities in equations (29e)–(29g). From the above discussions, the traceless feature of the multipole–multipole polarizabilities in equation (28) is the consequence of the traceless electric field gradients, attributed to the Laplace equation.

4. Numerical multipoles and (hyper)polarizabilities with finite field on ab initio calculations. A case study on carbon dioxide

Carbon dioxide has attracted much interest, not only for its importance in industry and environmental science, but also for its unique electronic properties as a quadrupole molecule. Buckingham et al [3, 7] had designed a four-wires condenser to measure the accurate quadrupole of CO₂, which laid the foundation for the electric-field-gradient-induced birefringence (EFGIB) experiments [17–19]. The mean and anisotropy dipole polarizability of CO₂ could be measured by the depolarization ratio of Rayleigh scattering [20]. Ritchie et al [21] utilized the Kerr effect to measure the mean second dipole hyperpolarizability, suggested by Buckingham and Pople [6]. On the other hand, higher order multipoles and (hyper)polarizabilities are difficult to directly measure from the experiments nowadays.

Computationally, Maroulis and co-workers have calculated the multipoles and (hyper)polarizabilities of CO₂ numerically against high-level ab initio calculations with both the charge perturbation method (CPM) [22] and the finite field method (FFM) [10, 23]. In the charge perturbation method, charge is placed at certain distance R from the molecule. Since all orders of electric field and gradients are generated by the charge, according to equation (14), the finite difference is based on the full Buckingham expansion. Moreover, for a polyatom molecule, the point multipole representations are valid only at R→∞. Thus, in the actual calculations, the multipoles and (hyper)polarizabilities are often calculated at several different R values, and then extrapolated to R→∞ value to refine the result [8]. Such difficulty of the CPM can be overcome by the FFM, in which the electric field or gradient tensors are manipulated directly to generate the desired uniform external field or gradients. In this study, we apply the FFM to calculate the multipoles and (hyper)polarizabilities of CO₂ numerically, as detailed in appendix B.

The equilibrium geometry of an electronic ground state CO₂, of D_∞h point group, has only one structural degree of freedom, i.e. the C–O bond length r_CO. In this study, we take r_CO=2.196 a₀ (Bohr), as deduced from precise experimental rotational constant of CO₂ [24]. The Cartesian coordinate system is set with the z-axis along the D_∞h symmetric axis and x, y, z-axes mutually perpendicular, corresponding to the principal axes of both α_αβ and Θ_αβ. In the principal frame, there are two independent components of α_αβ, one independent component of Θ_αβ, and three independent components of C_αβ,γδ, respectively [5]. With the above geometry and coordinate, the ab initio calculations, with finite external electric field or gradient applied, were performed at CCSD(T)/daug-cc-pvqz level of theory using the Gaussian16 package [15].

According to equation (6a), with U in equations (12) or (27), the dipole polarizability can be calculated numerically with finite difference, i.e. [25]:(30)$ \begin{eqnarray}\begin{array}{l}{\alpha }_{\alpha \beta }=\displaystyle \frac{4}{3}\times \displaystyle \frac{U(-{E}_{\alpha },{E}_{\beta })+U({E}_{\alpha },-{E}_{\beta })-U({E}_{\alpha },{E}_{\beta })-U(-{E}_{\alpha },-{E}_{\beta })}{4{E}_{\alpha }{E}_{\beta }}\\ \,-\,\displaystyle \frac{1}{3}\times \displaystyle \frac{U(-2{E}_{\alpha },2{E}_{\beta })+U(2{E}_{\alpha },-2{E}_{\beta })-U(2{E}_{\alpha },2{E}_{\beta })-U(-2{E}_{\alpha },-2{E}_{\beta })}{16{E}_{\alpha }{E}_{\beta }}+\,O({E}_{\alpha }^{4}),\end{array}\end{eqnarray}$in which U(E_α, E_β) denotes the electronic energy calculated at CCSD(T)/daug-cc-pvqz level of theory with finite field of E_α =E_β =±0.001 e/a₀² applied on α and β directions simultaneously. Equation (30) is derived via Taylor series with solely electric fields applied, as a special case of the Buckingham expansion. Since CO₂ is orientated along the principal axes of polarizability with the z-axis along the D_∞h symmetric axis, the calculations on α_zz and α_xx are simpler for the condition α=β. In this case, U(-E_α, E_α)=U₀ corresponds to the energy without field, and U(E_α, E_α)=U(2E_α), etc, so that five ab initio calculations are enough, instead of eight ab initio calculations, in order to calculate each of α_zz and α_xx with the FFM.

Similarly, the traced quadrupole and quadrupole–quadrupole polarizability can be calculated with:(31)$ \begin{eqnarray}\begin{array}{lll}{Q}_{\alpha \beta } & = & \displaystyle \frac{4}{3}\times \displaystyle \frac{U(-{E}_{\alpha \beta })-U({E}_{\alpha \beta })}{2{E}_{\alpha \beta }}\\ & & -\,\displaystyle \frac{1}{3}\times \displaystyle \frac{U(-2{E}_{\alpha \beta })-U(2{E}_{\alpha \beta })}{4{E}_{\alpha \beta }}+O({E}_{\alpha \beta }^{4})\end{array}\end{eqnarray}$(32)$ \begin{eqnarray}\begin{array}{l}{c}_{\alpha \beta ,\gamma \delta }=\displaystyle \frac{4}{3}\times \displaystyle \frac{U(-{E}_{\alpha \beta },{E}_{\gamma \delta })+U({E}_{\alpha \beta },-{E}_{\gamma \delta })-U({E}_{\alpha \beta },{E}_{\gamma \delta })-U(-{E}_{\alpha \beta },-{E}_{\gamma \delta })}{8{E}_{\alpha \beta }{E}_{\gamma \delta }}\\ \,-\,\displaystyle \frac{1}{3}\times \displaystyle \frac{U(-2{E}_{\alpha \beta },2{E}_{\gamma \delta })+U(2{E}_{\alpha \beta },-2{E}_{\gamma \delta })-U(2{E}_{\alpha \beta },2{E}_{\gamma \delta })-U(-2{E}_{\alpha \beta },-2{E}_{\gamma \delta })}{32{E}_{\alpha \beta }{E}_{\gamma \delta }}+\,O({E}_{\alpha \beta }^{4}).\end{array}\end{eqnarray}$

The detailed derivation of the above two finite field expressions are given by equations (B5) and (B8) in appendix A. The applied finite field gradient is E_αβ =0.0005 e/a₀³. Since the symmetric E_αβ tensor is manipulated in the Gaussian16 package [15] with E_βα =E_αβ =0.0005 e/a₀³ for α ≠ β and E_αα =0.001 e/a₀³ for α=β, the actual calculation with the FFM gives Q_αβ +Q_βα =2Q_αβ. Thus, the result needs to be further divided by two, as comparing equation (31) with equation (B5). For the same reason, the actual calculation with the FFM gives c_αβ,γδ +c_βα,γδ +c_αβ,δγ +c_βα,δγ =4c_αβ,γδ. Thus, the result needs to be further divided by four, as comparing equation (32) with equation (B8). The corresponding traceless quadrupole Θ_αβ and quadrupole–quadrupole polarizability C_αβ,γδ can be obtained via equations (22a) and (28c), respectively.

Tables 1–3 list the independent components of dipole polarizability, quadrupole and quadrupole–quadrupole polarizability, respectively. Comparisons with high-level ab initio calculations of both the CPM and the FFM, as well as the experimental measurements are also listed for dipole polarizability and quadrupoles in tables 1 and 2. It can be seen that the agreement of current study with previous computational and experimental studies is reasonable. Thus, the general Buckingham expansion in the traced form in equation (B1) is validated. The benefit of equation (B1) is that it can be applied to calculate all the multipoles and (hyper)polarizabilities with accuracy of O(Eⁿ) for arbitrary order n in a consistent manner with the FFM. Of course, the accuracy of the FFM is also highly dependent on the level of theory and basis set in the ab initio calculations, in order to get an accurate U(E). For further study, we plan to adopt high-level ab initio calculations and extrapolate to the complete basis set limit to get higher accuracy of U(E).


Table 1.
Table 1.The dipole polarizability of CO₂ in atomic units (a₀³).

	CCSD(T)a	MP4SDQb	CCSD(T)c	CCSD(T)d	Experiments
αzz	26.740	27.140	27.07	26.8184	27.2408e, 27.3056f, 27.2909g
αxxh	12.852	12.869	12.99	12.8966	13.0018e, 13.0162f, 13.0379g
$\bar{\alpha }$i	17.481	17.626	17.69	17.5372	17.7481e, 17.7793f, 17.7889g
Δαj	13.887	14.271	14.08	13.9217	14.2390e, 14.2894f, 14.2530g

^aThis work, the FFM against CCSD(T)/daug-cc-pvqz with 315 contracted Gaussian-type basis functions (CGTBFs).
^bReference [22], the CPM against MP4SDQ/[6s4p3d1f] with 120 CGTBFs.
^cReference [23], the FFM against CCSD(T)/[6s4p4d1f] with 135 CGTBFs.
^dReference [10], the FFM against CCSD(T)/[7s5p4d2f] with 168 CGTBFs.
^eReference [26], derived from the depolarization ratio of Rayleigh scattering.
^fReference [20], derived from the depolarization ratio of Rayleigh scattering.
^gReference [27], derived from the depolarization ratio of Rayleigh scattering.
^hFor CO₂ orientated with z-axis along the D_∞h symmetric axis, α_yy =α_xx.
ⁱIsotropic (average) dipole polarizability, $\bar{\alpha }$=(α_zz +2α_xx)/3.
^jAnisotropic dipole polarizability, Δα=α_zz−α_xx.

New window|CSV


Table 2.
Table 2.The traceless quadrupole of CO₂ in atomic units (e·a₀²).

	CCSD(T)a	MP4SDQb	CCSD(T)c	CCSD(T)d	Experiments
Θzzi	−3.163	−3.239	−3.19	−3.17	−3.0482e, −3.1873f, −3.1806g, −3.1895h

^aThis work.
^bReference [22].
^cReference [23].
^dReference [28], the FFM against CCSD(T)/[8s5p5d2f] with 186 CGTBFs.
^eReference [7], EFGIB experiment proposed by Buckingham.
^fReference [18], EFGIB experiment.
^gReference [19], EFGIB experiment.
^hReference [17], EFGIB experiment.
ⁱFor CO₂ orientated with z-axis along the D_∞h symmetric axis, the traceless condition in equation (25a) implies that Θ_xx =Θ_yy =−Θ_zz /2.

New window|CSV


Table 3.
Table 3.The traceless quadrupole–quadrupole polarizability of CO₂ in atomic units (a₀⁵).

	CCSD(T)a	MP4SDQb	CCSD(T)c	CCSD(T)d
Czz,zze	80.878	81.079	80.94	81.14
Cxx,xxf	34.082	33.021	34.13	33.97
Cxz,xzg	53.506	54.039	54.81	53.88
$\bar{C}$h	78.158	77.757	79.25	78.39

^aThis work.
^bReference [22].
^cReference [23].
^dReference [10].
^eFor CO₂ of D_∞h symmetry, the traceless condition in in equation (29c) implies that C_zz,xx =C_zz,yy =–C_zz,zz/2. Also, as discussed in section 2, the symmetry of subscripts implies that C_xx,zz =C_zz,xx =C_yy,zz =C_zz,yy.
^fC_xx,yy=C_yy,xx=–C_xx,xx−C_xx,zz.
^gC_yz,yz =C_xz,xz, C_xy,xy =C_xx,xx−C_zz,zz/4, and the values are invariant under the permutations of subscripts yz or xz or xy within the same group. All the other non-mentioned C_αβ,γδ’s are identically zeros.
^hIsotropic quadrupole–quadrupole polarizability $\bar{C}$=(C_zz,zz + 8C_xx,xx +8C_xz,xz)/10 [22].

New window|CSV

5. Summary

In this study, we derive the Buckingham expansion [3–5] in both traced and traceless forms. Specifically, with explicitly taking the Buckingham convention, aiming on the final form of the Buckingham expansion in its widely accepted form, we derive in section 2 the Buckingham expansion in the traced form using Taylor series. Based on the above derivation, a general Buckingham expansion in the traced form, equation (B1), is proposed. Based on equation (B1), numerical calculations of the multipoles and (hyper)polarizabilities with the finite field method, with accuracy of O(E⁴), in which E denotes the externally applied finite fields or gradients, is derived. The transformation from the traced multipoles and multipole–multipole polarizabilities to the corresponding traceless counterparts can be performed with an auxiliary traced electric field gradient, as described in section 3. The results of numerical calculations based on the finite field method against high-level ab initio calculations of CO₂, as a case study in section 4, are in good agreement with experimental results and previous theoretical calculations. We hope that the method in the current study can be applied to calculations of molecular/ionic multipoles and (hyper)polarizabilities, and further studies are in progress.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos. 21573112 and 21421001).

Conflict of interest

The authors declare that they have no conflict of interest.

Appendix A. Comparison with previous works

Buckingham proposed the Buckingham expansion in the traceless form [3–5], but did not provide a complete derivation to the best of our knowledge. McLean and Yoshimine (M&Y) [12], in communications with Buckingham, derived the Buckingham expansion in the traced form with multi-variable Taylor series, and then extended to the traceless form by solving the traceless multipoles using the Laplace equation, and solving the traceless multipole–multipole polarizabilities using the traceless features in equation (29). Our derivation also starts from Taylor series, as presented in section 2. On the other hand, the main difference is that the traceless multipoles and multipole–multipole polarizabilities are derived in section 3 in this work. Thus, equation (29), as the consequence of Laplace equation, is not needed in deriving the traceless multipole–multipole polarizabilities. Since M&Y’s derivation [12] was adopted in many studies in numerical calculations of molecular multipoles and (hyper)polarizabilities [23, 29], we give in this appendix mainly comparison with M&Y’s derivation [12].

Firstly, we need to clarify the convention in M&Y’s derivation [12]. For a charge distribution of ρ(r) centered at O, represented by point multipoles, the energy u, which is denoted as U in our study, is a function of electric potential, fields, and gradients, i.e. [12]:(A1)$ \begin{eqnarray}\begin{array}{lll}u & = & u\left(\varphi ,{E}_{x},{E}_{y},{E}_{z},{E}_{xx},{E}_{xy},{E}_{xz},{E}_{yy},{E}_{yz},{E}_{zz},\right.\\ & & {E}_{xxx},{E}_{xxy},{E}_{xxz},{E}_{xyy},{E}_{xyz},{E}_{xzz},{E}_{yyy},{E}_{yyz},\\ & & \left.{E}_{yzz},{E}_{zzz},{E}_{xxxx},\cdots \right),\end{array}\end{eqnarray}$in which the energy u₀ in zero external field is set to be zero without losing generality. It is notable that only half of the off-diagonal electric field gradient tensor components are included in the above equation. For example, E_xy is included in equation (A1), but E_yx is not included, since field gradient is symmetric and E_xy =E_yx. Thus, the convention in M&Y’s derivation [12] is that only one of the symmetric electric field gradient components is included in the energy expression. Similarly, by such convention, only one out of the three permutations of the xxy components is included for E_xxy, and one out of the six permutations of the xyz components is included for E_xyz, etc. Specifically, for an lth rank field gradient tensor, E_αβ⋯, only one out of the d_αβ⋯ permutations is included in equation (A1), with:(A2)$ \begin{eqnarray}{d}_{\alpha \beta \cdots }=\displaystyle \frac{l!}{{l}_{x}!{l}_{y}!{l}_{z}!},\end{eqnarray}$in which l_x, l_y, and l_z are the number of the repeated subscripts $\{x,y,z\},$ respectively, and l=l_x +l_y +l_z. By such convention, it is understood that only the independent field gradient components, i.e. 6 out of 9 E_αβ, 10 out of 27 E_αβγ, 15 out of 81 E_αβγδ, etc, are included in equation (A1). Generally, the number of independent components for an lth rank tensor is $\displaystyle {\sum }_{i=1}^{l+1}i.$ Such convention is in distinct difference than the Buckingham convention, in which all the permutations of E_αβ⋯ are included in the energy expression, as discussed in section 2.

Nevertheless, following the convention of independent field gradient components of M&Y [12], equation (A1) is then expanded via Taylor series in the space of multiple variables, i.e.:(A3)$ \begin{eqnarray}\begin{array}{lll}u & = & q\varphi +\displaystyle \displaystyle \sum _{\alpha }{u}_{\alpha }^{{\prime} }{E}_{\alpha }+\displaystyle \displaystyle \sum _{\alpha }\displaystyle \displaystyle \sum _{\beta \geqslant \alpha }{u}_{\alpha \beta }^{{\prime} }{E}_{\alpha \beta }\,+\,\displaystyle \displaystyle \sum _{\alpha }\displaystyle \displaystyle \sum _{\beta \geqslant \alpha }\displaystyle \displaystyle \sum _{\gamma \geqslant \beta }{u}_{\alpha \beta \gamma }^{{\prime} }{E}_{\alpha \beta \gamma }\\ & & +\,\displaystyle \displaystyle \sum _{\alpha }\displaystyle \displaystyle \sum _{\beta \geqslant \alpha }\displaystyle \displaystyle \sum _{\gamma \geqslant \beta }\displaystyle \displaystyle \sum _{\delta \geqslant \gamma }{u}_{\alpha \beta \gamma \delta }^{{\prime} }{E}_{\alpha \beta \gamma \delta }+\,\cdots +\,\displaystyle \frac{1}{2}\displaystyle \displaystyle \sum _{\alpha }\displaystyle \displaystyle \sum _{\beta }{u}_{\alpha \beta }^{{\prime\prime} }{E}_{\alpha }{E}_{\beta }\\ & & +\,\displaystyle \frac{1}{2}\displaystyle \displaystyle \sum _{\gamma }\displaystyle \displaystyle \sum _{\alpha }\displaystyle \displaystyle \sum _{\beta \geqslant \alpha }\left({u}_{\gamma ,\alpha \beta }^{{\prime\prime} }{E}_{\gamma }{E}_{\alpha \beta }+{u}_{\alpha \beta ,\gamma }^{{\prime\prime} }{E}_{\alpha \beta }{E}_{\gamma }\right)\\ \\ & & +\,\displaystyle \frac{1}{2}\displaystyle \displaystyle \sum _{\alpha }\displaystyle \displaystyle \sum _{\beta \geqslant \alpha }\displaystyle \displaystyle \sum _{\gamma }\displaystyle \displaystyle \sum _{\delta \geqslant \gamma }{u}_{\alpha \beta ,\gamma \delta }^{{\prime\prime} }{E}_{\alpha \beta }{E}_{\gamma \delta }\,+\,\displaystyle \frac{1}{2}\displaystyle \displaystyle \sum _{\delta }\displaystyle \displaystyle \sum _{\alpha }\displaystyle \displaystyle \sum _{\beta \geqslant \alpha }\displaystyle \displaystyle \sum _{\gamma \geqslant \beta }\left({u}_{\delta ,\alpha \beta \gamma }^{{\prime\prime} }{E}_{\delta }{E}_{\alpha \beta \gamma }\right.\\ & & \left.+{u}_{\alpha \beta \gamma ,\delta }^{{\prime\prime} }{E}_{\alpha \beta \gamma }{E}_{\delta }\right)+\cdots +\,\displaystyle \frac{1}{3!}\displaystyle \displaystyle \sum _{\alpha }\displaystyle \displaystyle \sum _{\beta }\displaystyle \displaystyle \sum _{\gamma }{u}_{\alpha \beta \gamma }^{\prime\prime\prime }{E}_{\alpha }{E}_{\beta }{E}_{\gamma }\\ & & +\,\displaystyle \frac{1}{3!}\displaystyle \displaystyle \sum _{\alpha }\displaystyle \displaystyle \sum _{\beta \geqslant \alpha }\displaystyle \displaystyle \sum _{\gamma }\displaystyle \displaystyle \sum _{\delta }\left({u}_{\alpha \beta ,\gamma \delta }^{\prime\prime\prime }{E}_{\alpha \beta }{E}_{\gamma }{E}_{\delta }\right.\\ & & \left.+\,{u}_{\gamma ,\alpha \beta ,\delta }^{\prime\prime\prime }{E}_{\gamma }{E}_{\alpha \beta }{E}_{\delta }+{u}_{\gamma \delta ,\alpha \beta }^{\prime\prime\prime }{E}_{\gamma }{E}_{\delta }{E}_{\alpha \beta }\right)+\,\cdots \\ & & +\,\displaystyle \frac{1}{4!}\displaystyle \displaystyle \sum _{\alpha }\displaystyle \displaystyle \sum _{\beta }\displaystyle \displaystyle \sum _{\gamma }\displaystyle \displaystyle \sum _{\delta }{u}_{\alpha \beta \gamma \delta }^{{\prime} \,^{\prime} \,^{\prime} \,^{\prime} }{E}_{\alpha }{E}_{\beta }{E}_{\gamma }{E}_{\delta }+\,\cdots ,\end{array}\end{eqnarray}$in which the summation symbol ∑ is explicitly written, because M&Y’s convention imposes the constraints on the range of subscripts for the lth rank field gradient tensors for l≥2. Apart from that, the exchange of the group of the subscripts, for example, in the constraint summation over ${u}_{\gamma ,\alpha \beta }^{{\prime\prime} }{E}_{\alpha }{E}_{\alpha \beta }+{u}_{\alpha \beta ,\gamma }^{{\prime\prime} }{E}_{\alpha \beta }{E}_{\gamma }$ in the second line of equation (A3), occurs as a consequence of the Taylor series of multi-variables. It is also notable that due to the inherent exchangeable groups of the sub-indices, the summations over ${u}_{\gamma \delta ,\alpha \beta }^{{\prime\prime} }{E}_{\gamma \delta }{E}_{\alpha \beta },$ are already included in the summations over ${u}_{\alpha \beta ,\gamma \delta }^{{\prime\prime} }{E}_{\alpha \beta }{E}_{\gamma \delta }$. Also, due to the free permutations of sub-indices γ and δ in ${u}_{\alpha \beta ,\gamma \delta }^{\prime\prime\prime },$ the summations over ${u}_{\alpha \beta ,\delta \gamma }^{\prime\prime\prime }{E}_{\alpha \beta }{E}_{\delta }{E}_{\gamma },$ have already been included in the summations over ${u}_{\alpha \beta ,\gamma \delta }^{\prime\prime\prime }{E}_{\alpha \beta }{E}_{\gamma }{E}_{\delta },$ and similar for the other two terms involving ${u}_{\gamma ,\alpha \beta ,\delta }^{\prime\prime\prime }$ and ${u}_{\gamma \delta ,\alpha \beta }^{\prime\prime\prime }.$ Thus, all the permutations of the groups of subscripts are allowed, and such convention is also in distinct difference than the Buckingham convention discussed in section 2, and a conversion to the Buckingham convention is needed.

Equation (A3) is the key to understand M&Y’s derivation [12] to solve the multipoles and (hyper)polarizabilities. Briefly, the approach taken in M&Y’s derivation [12] is that they firstly released the constraint on the independent electric field gradient variables in equation (A1), and then added a linear combinations of the traceless electric field gradients, i.e. equation (17), into equation (A3). This manipulation does not alter u, but equation (A3) can be transformed to the Buckingham expansion in the traceless form, i.e. equation (27). Next, the traceless features of the multipole–multipole polarizabilities, i.e. equation (29), were invited to solve the coefficients in the linear combination. By equating the coefficients, the individual Cartesian components of the multipoles and (hyper)polarizabilities were expressed in the derivatives of u w.r.t. electric fields or gradients [12].

Here, we take a simple route to get the same results in M&Y’s derivation [12], by taking advantage of some of the derived expressions in section 2 and section 3. Since the multi-variable Taylor series in equation (A3) ought to give the same energy as the Buckingham expansion, by comparing equivalent terms of equation (A3) with the Buckingham expansion in the traced form in equation (12), and considering the symmetric permutations, we have:(A4a)$ \begin{eqnarray}{u}_{\alpha }^{{\prime} }=\displaystyle \frac{\partial u}{\partial {E}_{\alpha }}=-{\mu }_{\alpha }\end{eqnarray}$(A4b)$ \begin{eqnarray}{u}_{\alpha \beta }^{{\prime} }=\displaystyle \frac{\partial u}{\partial {E}_{\alpha \beta }}=-\displaystyle \frac{{d}_{\alpha \beta }}{2}{Q}_{\alpha \beta }\end{eqnarray}$(A4c)$ \begin{eqnarray}{u}_{\alpha \beta \gamma }^{{\prime} }=\displaystyle \frac{\partial u}{\partial {E}_{\alpha \beta \gamma }}=-\displaystyle \frac{{d}_{\alpha \beta \gamma }}{3!}{O}_{\alpha \beta \gamma }\end{eqnarray}$(A4d)$ \begin{eqnarray}{u}_{\alpha \beta \gamma \delta }^{{\prime} }=\displaystyle \frac{\partial u}{\partial {E}_{\alpha \beta \gamma \delta }}=-\displaystyle \frac{{d}_{\alpha \beta \gamma \delta }}{4!}{H}_{\alpha \beta \gamma \delta }\end{eqnarray}$(A4e)$ \begin{eqnarray}{u}_{\alpha \beta }^{{\prime\prime} }=\displaystyle \frac{{\partial }^{2}u}{\partial {E}_{\alpha }\partial {E}_{\beta }}=-{\alpha }_{\alpha \beta }\end{eqnarray}$(A4f)$ \begin{eqnarray}{u}_{\gamma ,\alpha \beta }^{{\prime\prime} }={u}_{\alpha \beta ,\gamma }^{{\prime\prime} }=\displaystyle \frac{{\partial }^{2}u}{\partial {E}_{\gamma }\partial {E}_{\alpha \beta }}=-\displaystyle \frac{{d}_{\alpha \beta }}{2}{a}_{\gamma ,\alpha \beta }\end{eqnarray}$(A4g)$ \begin{eqnarray}{u}_{\alpha \beta ,\gamma \delta }^{{\prime\prime} }=\displaystyle \frac{{\partial }^{2}u}{\partial {E}_{\alpha \beta }\partial {E}_{\gamma \delta }}=-\displaystyle \frac{{d}_{\alpha \beta }{d}_{\gamma \delta }}{2}{c}_{\alpha \beta ,\gamma \delta }\end{eqnarray}$(A4h)$ \begin{eqnarray}{u}_{\delta ,\alpha \beta \gamma }^{{\prime\prime} }={u}_{\alpha \beta \gamma ,\delta }^{{\prime\prime} }=\displaystyle \frac{{\partial }^{2}u}{\partial {E}_{\delta }\partial {E}_{\alpha \beta \gamma }}=-\displaystyle \frac{{d}_{\alpha \beta \gamma }}{6}{d}_{\delta ,\alpha \beta \gamma }\end{eqnarray}$(A4i)$ \begin{eqnarray}{u}_{\alpha \beta \gamma }^{\prime\prime\prime }=\displaystyle \frac{{\partial }^{3}u}{\partial {E}_{\alpha }\partial {E}_{\beta }\partial {E}_{\gamma }}=-{\beta }_{\alpha \beta \gamma }\end{eqnarray}$(A4j)$ \begin{eqnarray}{u}_{\gamma \delta ,\alpha \beta }^{\prime\prime\prime }={u}_{\gamma ,\alpha \beta ,\delta }^{\prime\prime\prime }={u}_{\alpha \beta ,\gamma \delta }^{\prime\prime\prime }=\displaystyle \frac{{\partial }^{3}u}{\partial {E}_{\gamma }\partial {E}_{\delta }\partial {E}_{\alpha \beta }}=-\displaystyle \frac{{d}_{\alpha \beta }}{2}{b}_{\gamma \delta ,\alpha \beta }\end{eqnarray}$(A4k)$ \begin{eqnarray}{u}_{\alpha \beta \gamma \delta }^{{\prime} \,^{\prime} \,^{\prime} \,^{\prime} }=\displaystyle \frac{{\partial }^{4}u}{\partial {E}_{\alpha }\partial {E}_{\beta }\partial {E}_{\gamma }\partial {E}_{\delta }}=-{\gamma }_{\alpha \beta \gamma \delta }.\end{eqnarray}$It is notable that equations (A4a), (A4e), (A4i) and (A4k) give the dipole, polarizability, and first and second polarizability, as in equations (5a), (6a), (8a) and (8b). For the others, they give the relations between the derivatives of u and the traced multipoles and multipole–multipole polarizabilities in equations (5b)–(5d), (6b)–(6d) and (8c). We can next convert the traced multipoles and multipole–multipole polarizabilities to the traceless counterparts using equations (22) and (28), and the results are:(A5a)$ \begin{eqnarray}{{\rm{\Theta }}}_{\alpha \beta }=-3\displaystyle \frac{{u}_{\alpha \beta }^{{\prime} }}{{d}_{\alpha \beta }}+{u}_{\mu \mu }^{{\prime} }{\delta }_{\alpha \beta }\end{eqnarray}$(A5b)$ \begin{eqnarray}{{\rm{\Omega }}}_{\alpha \beta \gamma }=-15\displaystyle \frac{{u}_{\alpha \beta \gamma }^{{\prime} }}{{d}_{\alpha \beta \gamma }}+3\left(\displaystyle \frac{{u}_{\alpha \mu \mu }^{{\prime} }}{{d}_{\alpha \mu \mu }}{\delta }_{\beta \gamma }+\displaystyle \frac{{u}_{\mu \beta \mu }^{{\prime} }}{{d}_{\mu \beta \mu }}{\delta }_{\alpha \gamma }+\displaystyle \frac{{u}_{\mu \mu \gamma }^{{\prime} }}{{d}_{\mu \mu \gamma }}{\delta }_{\alpha \beta }\right)\end{eqnarray}$(A5c)$ \begin{eqnarray}\begin{array}{lll}{{\rm{\Phi }}}_{\alpha \beta \gamma \delta } & = & -105\displaystyle \frac{{u}_{\alpha \beta \gamma \delta }^{{\prime} }}{{d}_{\alpha \beta \gamma \delta }}+15\left(\displaystyle \frac{{u}_{\alpha \beta \mu \mu }^{{\prime} }}{{d}_{\alpha \beta \mu \mu }}{\delta }_{\gamma \delta }+\displaystyle \frac{{u}_{\alpha \mu \gamma \mu }^{{\prime} }}{{d}_{\alpha \mu \gamma \mu }}{\delta }_{\beta \delta }\right.\\ & & \left.+\,\displaystyle \frac{{u}_{\alpha \mu \mu \delta }^{{\prime} }}{{d}_{\alpha \mu \mu \delta }}{\delta }_{\beta \gamma }+\displaystyle \frac{{u}_{\mu \beta \gamma \mu }^{{\prime} }}{{d}_{\mu \beta \gamma \mu }}{\delta }_{\alpha \delta }+\displaystyle \frac{{u}_{\mu \beta \mu \delta }^{{\prime} }}{{d}_{\mu \beta \mu \delta }}{\delta }_{\alpha \gamma }+\displaystyle \frac{{u}_{\mu \mu \gamma \delta }^{{\prime} }}{{d}_{\mu \mu \gamma \delta }}{\delta }_{\alpha \beta }\right)\\ & & -\,3\left(\displaystyle \frac{{u}_{\mu \mu \nu \nu }^{{\prime} }}{{d}_{\mu \mu \nu \nu }}{\delta }_{\alpha \beta }{\delta }_{\gamma \delta }+\displaystyle \frac{{u}_{\mu \nu \mu \nu }^{{\prime} }}{{d}_{\mu \nu \mu \nu }}{\delta }_{\alpha \gamma }{\delta }_{\beta \delta }+\displaystyle \frac{{u}_{\mu \nu \nu \mu }^{{\prime} }}{{d}_{\mu \nu \nu \mu }}{\delta }_{\alpha \delta }{\delta }_{\beta \gamma }\right)\end{array}\end{eqnarray}$(A5d)$ \begin{eqnarray}{A}_{\gamma ,\alpha \beta }=-3\displaystyle \frac{{u}_{\gamma ,\alpha \beta }^{{\prime\prime} }}{{d}_{\alpha \beta }}+{u}_{\gamma ,\mu \mu }^{{\prime\prime} }{\delta }_{\alpha \beta }\end{eqnarray}$(A5e)$ \begin{eqnarray}{B}_{\gamma \delta ,\alpha \beta }=-3\displaystyle \frac{{u}_{\gamma \delta ,\alpha \beta }^{\prime\prime\prime }}{{d}_{\alpha \beta }}+{u}_{\gamma \delta ,\mu \mu }^{\prime\prime\prime }{\delta }_{\alpha \beta }\end{eqnarray}$(A5f)$ \begin{eqnarray}\begin{array}{lll}{C}_{\alpha \beta ,\gamma \delta } & = & -3\displaystyle \frac{{u}_{\alpha \beta ,\gamma \delta }^{{\prime\prime} }}{{d}_{\alpha \beta }{d}_{\gamma \delta }}+\displaystyle \frac{{u}_{\alpha \beta ,\mu \mu }^{{\prime\prime} }}{{d}_{\alpha \beta }}{\delta }_{\gamma \delta }+\displaystyle \frac{{u}_{\mu \mu ,\gamma \delta }^{{\prime\prime} }}{{d}_{\gamma \delta }}{\delta }_{\alpha \beta }\\ & & -\,\displaystyle \frac{1}{3}{u}_{\mu \mu ,\nu \nu }^{{\prime\prime} }{\delta }_{\alpha \beta }{\delta }_{\gamma \delta }\end{array}\end{eqnarray}$(A5g)$ \begin{eqnarray}\begin{array}{lll}{D}_{\delta ,\alpha \beta \gamma } & = & -15\displaystyle \frac{{u}_{\delta ,\alpha \beta \gamma }^{{\prime\prime} }}{{d}_{\alpha \beta \gamma }}+3\left(\displaystyle \frac{{u}_{\delta ,\alpha \mu \mu }^{{\prime\prime} }}{{d}_{\alpha \mu \mu }}{\delta }_{\beta \gamma }+\displaystyle \frac{{u}_{\delta ,\mu \beta \mu }^{{\prime\prime} }}{{d}_{\mu \beta \mu }}{\delta }_{\alpha \gamma }\right.\\ & & \left.+\,\displaystyle \frac{{u}_{\delta ,\mu \mu \gamma }^{{\prime\prime} }}{{d}_{\mu \mu \gamma }}{\delta }_{\alpha \beta }\right),\end{array}\end{eqnarray}$in which d_αβ is the permutation of the subscripts αβ. Note that $\tfrac{{u}_{\mu \mu }^{{\prime} }}{{d}_{\mu \mu }}=\tfrac{{u}_{xx}^{{\prime} }}{{d}_{xx}}+\tfrac{{u}_{yy}^{{\prime} }}{{d}_{yy}}+\tfrac{{u}_{xx}^{{\prime} }}{{d}_{xx}}={u}_{\mu \mu }^{{\prime} }$ in the last expression of equation (A5a) since d_xx =d_yy =d_zz =1, and similar for the last expressions of equations (A5d)–(A5f). Equation (A5) gives exactly the same traceless multipoles and multipole–multipole polarizabilities as in equation (20) of M&Y’s derivation [12], as can be easily checked. Moreover, equation (A5) is equivalent to equations (22) and (28), which give the multipoles and multipole–multipole polarizabilities in the Buckingham expansion in the traceless form, equation (27), and the permutation factor d_αβ⋯ in equation (A5) is given by equation (A2), in which only one of the symmetry d_αβ⋯ electric field gradient components is included. Of course, since equation (27) includes all the symmetric E_αβ⋯, the division of d_αβ⋯ is unavoidable for the conversion between equations (A3) and (12).

Finally, we note that some of the gaps in M&Y’s derivation [12] were filled by Applequist, who also adopted multi-variable Taylor series with all the electric field and gradient components included. The resultant expression is similar to equation (A3), but with all the Greek subscripts freely varying on $\{x,y,z\}.$ Thus, Applequist’s convention, with free permutations of subscripts and free exchanges of groups of subscripts, is also different from the Buckingham convention, and a conversion is needed. On the other hand, in the derivation we provide in section 2 in this study, the convention is explicitly designed at very beginning. Thus, the derivation is naturally driven to the final form of the Buckingham expansion in the traced form in equation (12), and no conversions are needed.

Appendix B. Finite field method based on the Buckingham expansion

Based on the discussions in section 2, we propose the following general expression of the Buckingham expansion in the traced form,(B1)$ \begin{eqnarray}U=q\varphi -\displaystyle \sum _{{n}_{\infty }=0}^{{N}_{\infty }}\cdots \displaystyle \sum _{{n}_{3}=0}^{{N}_{3}}\displaystyle \sum _{{n}_{2}=0}^{{N}_{2}}\displaystyle \sum _{{n}_{1}=0}^{{N}_{1}}{c}_{{n}_{1}{n}_{2}{n}_{3}\cdots }{P}_{{n}_{1}{n}_{2}{n}_{3}\cdots }{E}_{{n}_{1}{n}_{2}{n}_{3}\cdots },\end{eqnarray}$in which ${P}_{{n}_{1}{n}_{2}{n}_{3}\cdots }$ denotes either multipole or (hyper)polarizability in the traced form, with P₁=P_100⋯ =μ_α, P₀₁=P_0100⋯= Q_αβ, P₁₁=P_1100⋯=a_γ,αβ, P₀₂=P_0200⋯=c_αβ,γδ, etc, and P₀=P_000⋯=-U₀ as definition. It is understood that the exchange of the sub-indices, i.e. ${P}_{{n}_{2}{n}_{1}{n}_{3}\cdots },$ is not allowed. For example, for P₁₁ with n₁=1, and n₂=1, so that ${P}_{{n}_{1}{n}_{2}}$=a_γ,αβ is allowed, but ${P}_{{n}_{2}{n}_{1}}$=a_αβ,γ is not allowed. On the other hand, the permutations of the Greek subscripts in the same place of ${P}_{{n}_{1}{n}_{2}{n}_{3}\cdots }$ is allowed. For example, for P₀₂, both c_αβ,γδ and c_γδ,αβ are allowed because the groups of the Greek subscripts are inherently exchangeable, as discussed in section 2. Also, the permutation of subscripts in the same group is allowed. For example, a_γ,βα =a_γ,αβ are allowed. Such convention is of course consistent with the Buckingham convention, as discussed in section 2. In equation (B1):(B2)$ \begin{eqnarray}{c}_{{n}_{1}{n}_{2}{n}_{3}\cdots }=\displaystyle \frac{1}{{n}_{1}!}\displaystyle \frac{1}{2{!}^{{n}_{2}}}\displaystyle \frac{1}{3{!}^{{n}_{3}}}\cdots ,\,{E}_{{n}_{1}{n}_{2}{n}_{3}\cdots }={E}_{1}^{{n}_{1}}{E}_{2}^{{n}_{2}}{E}_{3}^{{n}_{3}}\cdots ,\end{eqnarray}$where ${E}_{l}^{{n}_{l}}$ denotes n_l power of the lth-rank electric field Cartesian tensor, with l Greek subscripts, so that ${E}_{1}^{1}={E}_{\alpha },$${E}_{1}^{2}={E}_{\alpha }{E}_{\beta },$${E}_{2}^{1}={E}_{\alpha \beta },$${E}_{2}^{2}={E}_{\alpha \beta }{E}_{\gamma \delta },$ and ${E}_{l}^{0}=1,$ etc.

Though equation (B1) can be expanded to N_∞, for practical usage only a small portion is enough. For example, let N₁=2, N₂=2, and N_l =0 for l>2, we have the truncated version of equation (B1), i.e.:
$ \begin{eqnarray*}\begin{array}{lll}U & = & q\varphi -\displaystyle \sum _{{n}_{2}=0}^{2}\displaystyle \sum _{{n}_{1}=0}^{2}{c}_{{n}_{1}{n}_{2}}{P}_{{n}_{1}{n}_{2}}{E}_{{n}_{1}{n}_{2}}\\ & = & q\varphi -{c}_{00}{P}_{00}{E}_{1}^{0}{E}_{2}^{0}-{c}_{10}{P}_{10}{E}_{1}^{1}{E}_{2}^{0}\\ & & -\,{c}_{20}{P}_{20}{E}_{1}^{2}{E}_{2}^{0}-{c}_{01}{P}_{01}{E}_{1}^{0}{E}_{2}^{1}-{c}_{11}{P}_{11}{E}_{1}^{1}{E}_{2}^{1}\\ & & -\,{c}_{21}{P}_{21}{E}_{1}^{2}{E}_{2}^{1}-{c}_{02}{P}_{02}{E}_{1}^{0}{E}_{2}^{2}-{c}_{12}{P}_{12}{E}_{1}^{1}{E}_{2}^{2}\\ & & -\,{c}_{22}{P}_{22}{E}_{1}^{2}{E}_{2}^{2}\\ & = & q\varphi +{U}_{0}-{\mu }_{\alpha }{E}_{\alpha }-\displaystyle \frac{1}{2}{\alpha }_{\alpha \beta }{E}_{\alpha }{E}_{\beta }-\displaystyle \frac{1}{2}{Q}_{\alpha \beta }{E}_{\alpha \beta }\\ & & -\,\displaystyle \frac{1}{2}{a}_{\gamma ,\alpha \beta }{E}_{\gamma }{E}_{\alpha \beta }\\ & & -\,\displaystyle \frac{1}{4}{b}_{\gamma \delta ,\alpha \beta }{E}_{\gamma }{E}_{\delta }{E}_{\alpha \beta }-\displaystyle \frac{1}{4}{c}_{\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }{E}_{\gamma \delta }\\ & & -\,\displaystyle \frac{1}{4}{c}_{\eta ,\alpha \beta ,\gamma \delta }^{(2)}{E}_{\eta }{E}_{\alpha \beta }{E}_{\gamma \delta }-\displaystyle \frac{1}{8}{c}_{\eta \xi ,\alpha \beta ,\gamma \delta }^{(3)}{E}_{\eta }{E}_{\xi }{E}_{\alpha \beta }{E}_{\gamma \delta },\end{array}\end{eqnarray*}$
which is comparable to equation (12), with an additional term ${P}_{22}={c}_{\eta \xi ,\alpha \beta ,\gamma \delta }^{(3)}$ called dipole–dipole–quadrupole–quadrupole polarizability. Note that the Einstein summation convention is applied on the repeated Greek subscripts in the above expression. For another example, when only the specific electric field ${E}_{{n}_{1}},$ i.e. (±E_α, ±E_β), is applied, the resultant Buckingham expansion is reduced to Taylor series of electric field, from which the finite difference expression of dipole polarizability α_αβ can be derived [25], as shown in equation (30) in section 4.

We next derive the finite field expression of the traced quadrupole Q_αβ using equation (B1). When only the specific electric field gradient of ${E}_{0{n}_{2}},$ i.e. ±E_αβ, is applied, only the corresponding ${P}_{0{n}_{2}}$ responses to the power of ${E}_{\alpha \beta }^{{n}_{2}},$ so that equation (B1) can be simplified to include the contribution solely from ±E_αβ, i.e.:(B3)$ \begin{eqnarray}\begin{array}{lll}U(\pm {E}_{\alpha \beta }) & = & -\displaystyle \sum _{{n}_{2}=0}^{\infty }{c}_{0{n}_{2}}{P}_{0{n}_{2}}{E}_{0{n}_{2}}\\ & = & {U}_{0}\mp \displaystyle \frac{1}{2}{Q}_{\alpha \beta }{E}_{\alpha \beta }-\displaystyle \frac{1}{4}{c}_{\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{2}\\ & & \mp \,\displaystyle \frac{1}{8}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{3}-\displaystyle \frac{1}{16}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{4}\\ & & +\,O({E}_{2}^{5}),\end{array}\end{eqnarray}$in which p_{αβ,αβ,αβ} is quadrupole–quadrupole–quadrupole polarizability, and p_{αβ,αβ,αβ,αβ} is quadrupole–quadrupole–quadrupole–quadrupole polarizability, respectively. Note that the Einstein summation convention is not applied to equation (B3) because the Cartesian component αβ is not dummy as it appears on the LHS of the expression. From equation (B3), Q_αβ can be obtained by two-side finite difference against E_αβ, i.e.:(B4a)$ \begin{eqnarray}\displaystyle \frac{U(-{E}_{\alpha \beta })-U({E}_{\alpha \beta })}{{E}_{\alpha \beta }}={Q}_{\alpha \beta }+\displaystyle \frac{1}{4}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{2}+O({E}_{2}^{4}),\end{eqnarray}$which gives Q_αβ with accuracy up to $O({E}_{2}^{2}).$ Doubling the electric field gradient, i.e. ±2E_αβ, we have:(B4b)$ \begin{eqnarray}\displaystyle \frac{U(-2{E}_{\alpha \beta })-U(2{E}_{\alpha \beta })}{2{E}_{\alpha \beta }}={Q}_{\alpha \beta }+{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{2}+O({E}_{2}^{4}).\end{eqnarray}$Combining equations (B4a) and (B4b) to eliminate the ${p}_{\alpha \beta ,\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{2}$ term, we have:(B5)$ \begin{eqnarray}\begin{array}{lll}{Q}_{\alpha \beta } & = & \displaystyle \frac{4}{3}\times \displaystyle \frac{U(-{E}_{\alpha \beta })-U({E}_{\alpha \beta })}{{E}_{\alpha \beta }}\\ & & -\,\displaystyle \frac{1}{3}\times \displaystyle \frac{U(-2{E}_{\alpha \beta })-U(2{E}_{\alpha \beta })}{2{E}_{\alpha \beta }}+O({E}_{2}^{4}).\end{array}\end{eqnarray}$The above expression is four-side finite difference against E_αβ with accuracy up to $O({E}_{2}^{4}),$ and is adopted in section 4 to calculate the traced quadrupole numerically.

Similarly, the quadrupole–quadrupole polarizability c_αβ,γδ can be obtained by applying two specific electric field gradients, (±E_αβ, ±E_γδ), simultaneously, and using equation (B1):(B6a)$ \begin{eqnarray}\begin{array}{lll}U(\pm {E}_{\alpha \beta },{E}_{\gamma \delta }) & = & -\displaystyle \sum _{{n}_{2}=0}^{\infty }{c}_{0{n}_{2}}{P}_{0{n}_{2}}{E}_{0{n}_{2}}\\ & = & {U}_{0}\mp \displaystyle \frac{1}{2}{Q}_{\alpha \beta }{E}_{\alpha \beta }-\displaystyle \frac{1}{2}{Q}_{\gamma \delta }{E}_{\gamma \delta }\\ & & -\,\displaystyle \frac{1}{4}{c}_{\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{2}\mp \displaystyle \frac{2}{4}{c}_{\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }{E}_{\gamma \delta }\\ & & -\,\displaystyle \frac{1}{4}{c}_{\gamma \delta ,\gamma \delta }{E}_{\gamma \delta }^{2}\mp \displaystyle \frac{1}{8}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{3}\\ & & -\,\displaystyle \frac{3}{8}{p}_{\alpha \beta ,\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }^{2}{E}_{\gamma \delta }\\ & & \mp \,\displaystyle \frac{3}{8}{p}_{\alpha \beta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }{E}_{\gamma \delta }^{2}\\ & & -\,\displaystyle \frac{1}{8}{p}_{\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\gamma \delta }^{3}-\displaystyle \frac{1}{16}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{4}\\ & & \mp \,\displaystyle \frac{4}{16}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }^{3}{E}_{\gamma \delta }\\ & & -\,\displaystyle \frac{6}{16}{p}_{\alpha \beta ,\alpha \beta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }^{2}{E}_{\gamma \delta }^{2}\\ & & \mp \,\displaystyle \frac{4}{16}{p}_{\alpha \beta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }{E}_{\gamma \delta }^{3}\\ & & -\,\displaystyle \frac{1}{16}{p}_{\gamma \delta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\gamma \delta }^{4}\\ & & \mp \,\displaystyle \frac{1}{32}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{5}\\ & & -\,\displaystyle \frac{5}{32}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }^{4}{E}_{\gamma \delta }\\ & & \mp \,\displaystyle \frac{10}{32}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }^{3}{E}_{\gamma \delta }^{2}\\ & & -\,\displaystyle \frac{10}{32}{p}_{\alpha \beta ,\alpha \beta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }^{2}{E}_{\gamma \delta }^{3}\\ & & \mp \,\displaystyle \frac{5}{32}{p}_{\alpha \beta ,\gamma \delta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }{E}_{\gamma \delta }^{4}\\ & & -\,\displaystyle \frac{1}{32}{p}_{\gamma \delta ,\gamma \delta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\gamma \delta }^{5}+O({E}_{2}^{6})\end{array}\end{eqnarray}$(B6b)$ \begin{eqnarray}\begin{array}{lll}U(\pm {E}_{\alpha \beta },-{E}_{\gamma \delta }) & = & -\displaystyle \sum _{{n}_{2}=0}^{\infty }{c}_{0{n}_{2}}{P}_{0{n}_{2}}{E}_{0{n}_{2}}\\ & = & {U}_{0}\mp \displaystyle \frac{1}{2}{Q}_{\alpha \beta }{E}_{\alpha \beta }+\displaystyle \frac{1}{2}{Q}_{\gamma \delta }{E}_{\gamma \delta }\\ & & -\,\displaystyle \frac{1}{4}{c}_{\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{2}\pm \displaystyle \frac{2}{4}{c}_{\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }{E}_{\gamma \delta }\\ & & -\,\displaystyle \frac{1}{4}{c}_{\gamma \delta ,\gamma \delta }{E}_{\gamma \delta }^{2}\mp \displaystyle \frac{1}{8}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{3}\\ & & +\,\displaystyle \frac{3}{8}{p}_{\alpha \beta ,\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }^{2}{E}_{\gamma \delta }\\ & & \mp \,\displaystyle \frac{3}{8}{p}_{\alpha \beta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }{E}_{\gamma \delta }^{2}+\displaystyle \frac{1}{8}{p}_{\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\gamma \delta }^{3}\\ & & -\,\displaystyle \frac{1}{16}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{4}\\ & & \pm \,\displaystyle \frac{4}{16}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }^{3}{E}_{\gamma \delta }\\ & & -\,\displaystyle \frac{6}{16}{p}_{\alpha \beta ,\alpha \beta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }^{2}{E}_{\gamma \delta }^{2}\\ & & \pm \,\displaystyle \frac{4}{16}{p}_{\alpha \beta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }{E}_{\gamma \delta }^{3}\\ & & -\,\displaystyle \frac{1}{16}{p}_{\gamma \delta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\gamma \delta }^{4}\\ & & \mp \,\displaystyle \frac{1}{32}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\alpha \beta ,\alpha \beta }{E}_{\alpha \beta }^{5}\\ & & +\,\displaystyle \frac{5}{32}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }^{4}{E}_{\gamma \delta }\\ & & \mp \,\displaystyle \frac{10}{32}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }^{3}{E}_{\gamma \delta }^{2}\\ & & +\,\displaystyle \frac{10}{32}{p}_{\alpha \beta ,\alpha \beta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }^{2}{E}_{\gamma \delta }^{3}\\ & & \mp \,\displaystyle \frac{5}{32}{p}_{\alpha \beta ,\gamma \delta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\alpha \beta }{E}_{\gamma \delta }^{4}\\ & & +\,\displaystyle \frac{1}{32}{p}_{\gamma \delta ,\gamma \delta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\gamma \delta }^{5}+O({E}_{2}^{6}),\end{array}\end{eqnarray}$in which the coefficients in the denominator are determined according to equation (B2), while the coefficients in the numerator are determined according to permutation of the groups of the subscripts, which are inherently exchangeable, as discussed in section 2. From equation (B6), c_αβ,γδ can be obtained by four-side finite difference against E_αβ and E_γδ, i.e.:(B7a)$ \begin{eqnarray}\begin{array}{l}\displaystyle \frac{U(-{E}_{\alpha \beta },{E}_{\gamma \delta })+U({E}_{\alpha \beta },-{E}_{\gamma \delta })-U({E}_{\alpha \beta },{E}_{\gamma \delta })-U(-{E}_{\alpha \beta },-{E}_{\gamma \delta })}{2{E}_{\alpha \beta }{E}_{\gamma \delta }}\\ \,=\,{c}_{\alpha \beta ,\gamma \delta }+\displaystyle \frac{8}{16}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }^{2}+\displaystyle \frac{8}{16}{p}_{\alpha \beta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\gamma \delta }^{2}+O({E}_{2}^{4}),\end{array}\end{eqnarray}$which gives c_αβ,γδ with accuracy up to $O({E}_{2}^{2}).$ Doubling the electric field gradient, i.e. ±2E_αβ and ±2E_γδ, we have:(B7b)$ \begin{eqnarray}\begin{array}{l}\displaystyle \frac{U(-2{E}_{\alpha \beta },2{E}_{\gamma \delta })+U(2{E}_{\alpha \beta },-2{E}_{\gamma \delta })-U(2{E}_{\alpha \beta },2{E}_{\gamma \delta })-U(-2{E}_{\alpha \beta },-2{E}_{\gamma \delta })}{8{E}_{\alpha \beta }{E}_{\gamma \delta }}\\ \,=\,{c}_{\alpha \beta ,\gamma \delta }+\displaystyle \frac{32}{16}{p}_{\alpha \beta ,\alpha \beta ,\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }^{2}+\displaystyle \frac{32}{16}{p}_{\alpha \beta ,\gamma \delta ,\gamma \delta ,\gamma \delta }{E}_{\gamma \delta }^{2}+O({E}_{2}^{4}).\end{array}\end{eqnarray}$Combining equations (B7a) and (B7b) to eliminate the ${E}_{2}^{2}$ terms, we have:(B8)$ \begin{eqnarray}\begin{array}{l}{c}_{\alpha \beta ,\gamma \delta }=\displaystyle \frac{4}{3}\times \displaystyle \frac{U(-{E}_{\alpha \beta },{E}_{\gamma \delta })+U({E}_{\alpha \beta },-{E}_{\gamma \delta })-U({E}_{\alpha \beta },{E}_{\gamma \delta })-U(-{E}_{\alpha \beta },-{E}_{\gamma \delta })}{2{E}_{\alpha \beta }{E}_{\gamma \delta }}\\ \,-\,\displaystyle \frac{1}{3}\times \displaystyle \frac{U(-2{E}_{\alpha \beta },2{E}_{\gamma \delta })+U(2{E}_{\alpha \beta },-2{E}_{\gamma \delta })-U(2{E}_{\alpha \beta },2{E}_{\gamma \delta })-U(-2{E}_{\alpha \beta },-2{E}_{\gamma \delta })}{8{E}_{\alpha \beta }{E}_{\gamma \delta }}\\ \,+\,O({E}_{2}^{4}).\end{array}\end{eqnarray}$The above expression is eight-side finite difference against E_αβ and E_γδ with accuracy up to $O({E}_{2}^{4}),$ and is adopted in section 4 to calculate the traced quadrupole–quadrupole polarizability numerically. Note that the actual calculations for certain component may be simplified because the simultaneously applied two electric field gradients in opposite directions cancel each other. For example, using equation (B8), c_zz,zz can be calculated by:
$ \begin{eqnarray*}\begin{array}{lll}{c}_{zz,zz} & = & \displaystyle \frac{4}{3}\times \displaystyle \frac{2{U}_{0}-U(2{E}_{zz})-U(-2{E}_{zz})}{2{E}_{zz}^{2}}\\ & & -\,\displaystyle \frac{1}{3}\times \displaystyle \frac{2{U}_{0}-U(4{E}_{zz})-U(-4{E}_{zz})}{8{E}_{zz}^{2}}+O({E}_{zz}^{4}),\end{array}\end{eqnarray*}$
in which U₀ is the energy with zero external gradient because the simultaneously applied (+E_zz, −E_zz) cancel each other. Therefore, only five, instead of eight, ab initio calculations are needed to calculate c_zz,zz numerically with the finite field method.

The above derivation on Q_αβ and c_αβ,γδ are general and can be easily applied to calculate other traced multipoles or (hyper)polarizabilities of interest in the frame of the finite field method. To obtain accurate polarizabilities/hyperpolarizabilities, wavefunction-based methods such as CCSD(T) using the finite field approach are necessary. And for large π-conjugated systems such as polydiacetylene (PDA) and polybutatriene (PBT) as shown in previous works [30, 31], DFT methods with improved exchange-correlation (XC) functionals are also good choices when considering the balance between accuracy and computational cost. It should be noted that our finite field method derived in appendix A could be applied in conjunction with both wavefunction-based and DFT methods to obtain higher accuracy. The transformations to the corresponding traceless multipoles or (hyper)polarizabilities are discussed in section 3 as well as appendix B.

Appendix C. General Buckingham expansion in the traceless form

The general formula of the electric field gradients is [13, 32]:(C1)$ \begin{eqnarray}\begin{array}{lll}{E}_{M} & = & {E}_{{\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{M}}\\ & = & e{R}^{-\left(2M+1\right)}\displaystyle \sum _{L=0}^{\left[\displaystyle \frac{M}{2}\right]}{\left(-1\right)}^{L}\left[2\left(M-L\right)-1\right]!!{R}^{2L}{\left({R}^{M-2L}{\delta }^{L}\right)}_{{\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{M}}\\ & = & \displaystyle \sum _{L=0}^{\left[\displaystyle \frac{M}{2}\right]}{\left(-1\right)}^{L}\displaystyle \frac{\left[2\left(M-L\right)-1\right]!!}{\left(2M-1\right)!!}{\left(F{\delta }^{L}\right)}_{{\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{M}},\end{array}\end{eqnarray}$in which M stands for the rank of the traceless electric field gradients E_M =${E}_{{\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{M}},$F=F_M =${F}_{{\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{M}}$ denotes the corresponding auxiliary traced electric field gradients, [M/2] represents the largest integer not exceeding M/2. And in equation (C1):(C2)$ \begin{eqnarray}\begin{array}{lll}{\left({R}^{M-2L}{\delta }^{L}\right)}_{{\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{M}} & = & \sum {\delta }_{{\alpha }_{1}{\alpha }_{2}}{\delta }_{{\alpha }_{3}{\alpha }_{4}}\cdots {\delta }_{{\alpha }_{2L-1}{\alpha }_{2L}}{R}_{{\alpha }_{2L+1}}{R}_{{\alpha }_{2L+2}}\cdots {R}_{{\alpha }_{M-1}}{R}_{{\alpha }_{M}}\\ {\left(F{\delta }^{L}\right)}_{{\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{M}} & = & \sum {\delta }_{{\alpha }_{1}{\alpha }_{2}}{\delta }_{{\alpha }_{3}{\alpha }_{4}}\cdots {\delta }_{{\alpha }_{2L-1}{\alpha }_{2L}}{F}_{\mu \mu \nu \nu \cdots \xi \xi {\alpha }_{2L+1}{\alpha }_{2L+2}\cdots {\alpha }_{M-1}{\alpha }_{M}},\end{array}\end{eqnarray}$in which the sum is over all the different permutations of the subscripts α₁, α₂, ··· , α_M, including N_M,L =M!/[2^LL!(M − 2L)!] terms.

Based on equations (C1) and (C2), the general formula of equations (21a)–(21c) is given as:(C3)$ \begin{eqnarray}\begin{array}{lll}{J}_{{E}_{M},{F}_{M}} & = & {J}_{{\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{M},{\beta }_{1}{\beta }_{2}\cdots {\beta }_{M}}\\ & = & \displaystyle \frac{\partial {E}_{{\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{M}}}{\partial {F}_{{\beta }_{1}{\beta }_{2}\cdots {\beta }_{M}}}\\ & = & \displaystyle \sum _{L=0}^{\left[\displaystyle \frac{M}{2}\right]}{\left(-1\right)}^{L}\displaystyle \frac{\left[2\left(M-L\right)-1\right]!!}{\left(2M-1\right)!!}{\left({\delta }_{{\alpha }_{o}{\beta }_{o}}^{M-2L}{\delta }_{{\alpha }_{p}{\alpha }_{q}}^{L}{\delta }_{{\beta }_{p}{\beta }_{q}}^{L}\right)}_{\begin{array}{l}{\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{M}\\ {\beta }_{1}{\beta }_{2}\cdots {\beta }_{M}\end{array}}\end{array}\end{eqnarray}$and(C4)$ \begin{eqnarray}\begin{array}{lll}{\left({\delta }_{{\alpha }_{o}{\beta }_{o}}^{M-2L}{\delta }_{{\alpha }_{p}{\alpha }_{q}}^{L}{\delta }_{{\beta }_{p}{\beta }_{q}}^{L}\right)}_{\begin{array}{l}{\alpha }_{1}{\alpha }_{2}\cdots {\alpha }_{M}\\ {\beta }_{1}{\beta }_{2}\cdots {\beta }_{M}\end{array}} & = & \sum {\delta }_{{\alpha }_{1}{\alpha }_{2}}{\delta }_{{\beta }_{1}{\beta }_{2}}{\delta }_{{\alpha }_{3}{\alpha }_{4}}{\delta }_{{\beta }_{3}{\beta }_{4}}\cdots \\ & & \times \,{\delta }_{{\alpha }_{2L-1}{\alpha }_{2L}}{\delta }_{{\beta }_{2L-1}{\beta }_{2L}}{\delta }_{{\alpha }_{2L+1}{\beta }_{2L+1}}\\ & & \times \,{\delta }_{{\alpha }_{2L+2}{\beta }_{2L+2}}\cdots {\delta }_{{\alpha }_{M}{\beta }_{M}},\end{array}\end{eqnarray}$in which the sum is over the all the different permutations of the subscripts α₁, α₂, ··· , α_M, and the permutations of the subscripts β₁, β₂, ··· , β_M need to follow the same order.

With the help of equations (C3) and (C4), the general formula from traced to traceless multipoles and (hyper)polarizabilities is given as:(C5)$ \begin{eqnarray}{{\rm{\Xi }}}_{{n}_{1}{n}_{2}{n}_{3}\cdots }={\Im }_{{n}_{1}{n}_{2}{n}_{3}\cdots }{P}_{{n}_{1}{n}_{2}{n}_{3}\cdots },\end{eqnarray}$in which ${P}_{{n}_{1}{n}_{2}{n}_{3}\cdots }$ is the traced multipole or (hyper)polarizability defined in the equation (B1), and ${{\rm{\Xi }}}_{{n}_{1}{n}_{2}{n}_{3}\cdots }$ denotes the corresponding traceless multipole or (hyper)polarizability, with Ξ₀₁=_0100⋯=Θ_αβ, Ξ₁₁=Ξ_1100⋯=A_γ,αβ, Ξ₀₂= Ξ_0200⋯=C_αβ,γδ, etc, while Ξ₁=Ξ_100⋯= P_100⋯=μ_α, Ξ₂=Ξ_200⋯= P_200⋯=α_αβ, Ξ₃=Ξ_300⋯= P_300⋯=β_αβγ, Ξ₄=Ξ_400⋯= P_400⋯=γ_αβγδ, and Ξ₀=Ξ_000⋯= P_000⋯= −U₀ as definition. In equation (C5), the operator ${\Im }_{{n}_{1}{n}_{2}{n}_{3}\cdots }$ is:(C6)$ \begin{eqnarray}\begin{array}{lll}{\Im }_{{n}_{1}{n}_{2}{n}_{3}\cdots } & = & \underset{M}{\overset{\left\{M| {n}_{M}\ne 0,M=2,3,\cdots \right\}}{{\rm{\Pi }}}}{f}_{M}{\left({J}_{{E}_{M},{F}_{M}}\right)}^{{n}_{M}}\\ & = & \underset{M}{\overset{\left\{M| {n}_{M}\ne 0,M=2,3,\cdots \right\}}{{\rm{\Pi }}}}\displaystyle \frac{\left(2M-1\right)!!}{M!}{\left({J}_{{E}_{M},{F}_{M}}\right)}^{{n}_{M}},\end{array}\end{eqnarray}$in which just one corresponding f_M factor in the equation (23) is needed in the operator when n_M ≠ 0, while ${J}_{{E}_{M},{F}_{M}}$ need to be carried out n_M times to satisfy the traceless conditions.

Based on the discussions in section 3 and equations (C1)–(C6), we propose the following general expression of the Buckingham expansion in the traceless form:(C7)$ \begin{eqnarray}U=q\varphi -\displaystyle \sum _{{n}_{\infty }=0}^{{N}_{\infty }}\cdots \displaystyle \sum _{{n}_{3}=0}^{{N}_{3}}\displaystyle \sum _{{n}_{2}=0}^{{N}_{2}}\displaystyle \sum _{{n}_{1}=0}^{{N}_{1}}{f}_{{n}_{1}{n}_{2}{n}_{3}\cdots }{c}_{{n}_{1}{n}_{2}{n}_{3}\cdots }{{\rm{\Xi }}}_{{n}_{1}{n}_{2}{n}_{3}\cdots }{E}_{{n}_{1}{n}_{2}{n}_{3}\cdots },\end{eqnarray}$in which ${c}_{{n}_{1}{n}_{2}{n}_{3}\cdots }$ and ${E}_{{n}_{1}{n}_{2}{n}_{3}\cdots }$ have been defined in equation (B2). And in equation (C7):(C8)$ \begin{eqnarray}\begin{array}{lll}{f}_{{n}_{1}{n}_{2}{n}_{3}\cdots } & = & \underset{M}{\overset{\left\{M| {n}_{M}\ne 0,M=2,3,\cdots \right\}}{{\rm{\Pi }}}}\displaystyle \frac{M!}{(2M-1)!!}\\ & = & \underset{M}{\overset{\left\{M| {n}_{M}\ne 0,M=2,3,\cdots \right\}}{{\rm{\Pi }}}}\displaystyle \frac{1}{{f}_{M}},\end{array}\end{eqnarray}$where ${f}_{{n}_{1}{n}_{2}{n}_{3}\cdots }$ is the conversion coefficient from traced to traceless multipole or multipole–multipole polarizability tensor, defining that when n₂ =n₃ =···=0, ${f}_{{n}_{1}{n}_{2}{n}_{3}\cdots }=1.$

It is understood that in equation (C7), the exchange of the sub-indices, i.e. ${{\rm{\Xi }}}_{{n}_{2}{n}_{1}{n}_{3}\cdots },$ is not allowed. For example, for Ξ₁₁ with n₁=1, and n₂=1, so that ${{\rm{\Xi }}}_{{n}_{1}{n}_{2}}$=A_γ,αβ is allowed, but ${{\rm{\Xi }}}_{{n}_{2}{n}_{1}}$=A_αβ,γ is not allowed. On the other hand, the permutations of the Greek subscripts in the same place of ${{\rm{\Xi }}}_{{n}_{1}{n}_{2}{n}_{3}\cdots }$ are allowed. For example, for Ξ₀₂, both C_αβ,γδ and C_γδ,αβ are allowed because the groups of the Greek subscripts are inherently exchangeable, as discussed in the text. Also, the permutation of subscripts in the same group is allowed. For example, A_γ,βα =A_γ,αβ are allowed. Such convention is consistent with the Buckingham convention, as discussed in the text. It is noted that the equation (C8) means that when n_M ≠ 0, just one corresponding 1/f_M factor is needed in the transform, as shown in equations (26a)–(26c). For example, when n₁=0, n₂=2, the corresponding energy term transform in equation (26a) from traced to traceless expression is $-\tfrac{1}{4}{c}_{\alpha \beta ,\gamma \delta }{E}_{\gamma \delta }{E}_{\alpha \beta }$$=-\tfrac{1}{{f}_{{\rm{\Theta }}}}\tfrac{1}{4}{f}_{{\rm{\Theta }}}{c}_{\alpha \beta ,\gamma \delta }{E}_{\gamma \delta }{E}_{\alpha \beta }$$=-\tfrac{1}{6}{C}_{\alpha \beta ,\gamma \delta }{E}_{\gamma \delta }{E}_{\alpha \beta },$ in which f_Θ appears just once.

Though equation (C7) can be expanded to N_∞, for practice only a small portion is enough. For example, let N₁=2, N₂=2, and N_l =0 for l>2, we have the truncated version of equation (C7), i.e.:
$ \begin{eqnarray*}\begin{array}{lll}U & = & q\varphi -\displaystyle \sum _{{n}_{2}=0}^{2}\displaystyle \sum _{{n}_{1}=0}^{2}{f}_{{n}_{1}{n}_{2}}{c}_{{n}_{1}{n}_{2}}{{\rm{\Xi }}}_{{n}_{1}{n}_{2}}{E}_{{n}_{1}{n}_{2}}\\ & = & q\varphi -{f}_{00}{c}_{00}{{\rm{\Xi }}}_{00}{E}_{1}^{0}{E}_{2}^{0}-{f}_{10}{c}_{10}{{\rm{\Xi }}}_{10}{E}_{1}^{1}{E}_{2}^{0}\\ & & -\,{f}_{20}{c}_{20}{{\rm{\Xi }}}_{20}{E}_{1}^{2}{E}_{2}^{0}-{f}_{01}{c}_{01}{{\rm{\Xi }}}_{01}{E}_{1}^{0}{E}_{2}^{1}\\ & & -\,{f}_{11}{c}_{11}{{\rm{\Xi }}}_{11}{E}_{1}^{1}{E}_{2}^{1}-{f}_{21}{c}_{21}{{\rm{\Xi }}}_{21}{E}_{1}^{2}{E}_{2}^{1}\\ & & -\,{f}_{02}{c}_{02}{{\rm{\Xi }}}_{02}{E}_{1}^{0}{E}_{2}^{2}-{f}_{12}{c}_{12}{{\rm{\Xi }}}_{12}{E}_{1}^{1}{E}_{2}^{2}\\ & & -\,{f}_{22}{c}_{22}{{\rm{\Xi }}}_{22}{E}_{1}^{2}{E}_{2}^{2}\\ & = & q\varphi +{U}_{0}-{\mu }_{\alpha }{E}_{\alpha }-\displaystyle \frac{1}{2}{\alpha }_{\alpha \beta }{E}_{\alpha }{E}_{\beta }\\ & & -\,\displaystyle \frac{1}{3}{{\rm{\Theta }}}_{\alpha \beta }{E}_{\alpha \beta }-\displaystyle \frac{1}{3}{A}_{\gamma ,\alpha \beta }{E}_{\gamma }{E}_{\alpha \beta }\\ & & -\,\displaystyle \frac{1}{6}{B}_{\gamma \delta ,\alpha \beta }{E}_{\gamma }{E}_{\delta }{E}_{\alpha \beta }-\displaystyle \frac{1}{6}{C}_{\alpha \beta ,\gamma \delta }{E}_{\alpha \beta }{E}_{\gamma \delta }\\ & & -\,\displaystyle \frac{1}{6}{C}_{\eta ,\alpha \beta ,\gamma \delta }^{(2)}{E}_{\eta }{E}_{\alpha \beta }{E}_{\gamma \delta }\\ & & -\,\displaystyle \frac{1}{12}{C}_{\eta \xi ,\alpha \beta ,\gamma \delta }^{(3)}{E}_{\eta }{E}_{\xi }{E}_{\alpha \beta }{E}_{\gamma \delta },\end{array}\end{eqnarray*}$
which is consistent with the equation (27).

Other conventions may be adopted in the energy expansion besides the Buckingham convention discussed above. For example, Applequist [13] defined the general traced multipoles or (hyper)polarizabilities as:(C9)$ \begin{eqnarray}{{\boldsymbol{p}}}^{\left({k}_{1},{k}_{2},\cdots ,{k}_{t}\right)}=-\left(\displaystyle \frac{{\partial }^{t}U}{\partial {{\boldsymbol{E}}}^{\left({k}_{1}\right)}\partial {{\boldsymbol{E}}}^{\left({k}_{2}\right)}\cdots \partial {{\boldsymbol{E}}}^{\left({k}_{t}\right)}}\right),\end{eqnarray}$in which U is the total electrostatic energy, ${{\boldsymbol{E}}}^{\left({k}_{i}\right)}$ represents electric field or gradients of rank k_t which is the same as our traceless definition, i.e. equation (C1) and when t=1, ${{\boldsymbol{p}}}^{\left({k}_{1}\right)}$ is the traced multipole with p⁽¹⁾=μ, p⁽²⁾=1/2Q, p⁽³⁾= 1/6O, p⁽⁴⁾=1/24H, ···. Applequist gave his energy expansion as:(C10)$ \begin{eqnarray}U=q\varphi -\displaystyle \sum _{t=1}^{\infty }\displaystyle \frac{1}{t!}\displaystyle \displaystyle \sum _{{k}_{1},{k}_{2},\cdots ,{k}_{t}}{{\boldsymbol{p}}}^{\left({k}_{1},{k}_{2},\cdots ,{k}_{t}\right)}{{\boldsymbol{E}}}^{\left({k}_{1}\right)}{{\boldsymbol{E}}}^{\left({k}_{2}\right)}\cdots {{\boldsymbol{E}}}^{\left({k}_{t}\right)},\end{eqnarray}$in which k₁, k₂, ···, k_t could traverse from 1 to ∞. The equation (C10) is just the multi-variable Taylor expansion where the variables are the electric field and gradients E^(k). The convention adopted by Applequist in equation (C10) is that the permutations of groups k₁, k₂, ···, k_t are permitted, which is different from the Buckingham convention, and the permutations of subscripts inside the group are also permitted which is the same as the Buckingham convention. For example, there exists both p^(1,2) and p^(2,1) terms in the equation (C10), while there just exists a_γ,αβ, not a_αβ,γ in the equation (B1). It should be noted that the energy in equations (B1) and (C10) are equal. By comparing equations (B1) and (C10), we could find the relationships between Applequist’s and our traced definition, i.e.:(C11)$ \begin{eqnarray}{{\bf{p}}}^{\left({k}_{1},{k}_{2},\cdots ,{k}_{t}\right)}=\displaystyle \prod _{M}^{\left\{M| {n}_{M}\ne 0,M=2,3,\cdots \right\}}\displaystyle \frac{{n}_{M}!}{M{!}^{{n}_{M}}}{P}_{{n}_{1}{n}_{2}{n}_{3}\cdots },\end{eqnarray}$in which ${P}_{{n}_{1}{n}_{2}{n}_{3}\cdots }$ is the traced multipole or (hyper)polarizability defined in the equation (B1) and n₁+n₂+ n₃+ ···=t. Then, Applequist defined the general traceless multipole and (hyper)polarizabilities as:(C12)$ \begin{eqnarray}{{\boldsymbol{\xi }}}^{\left({k}_{1},{k}_{2},\cdots ,{k}_{t}\right)}={{T}}_{{k}_{1}}{{T}}_{{k}_{2}}\cdots {{T}}_{{k}_{t}}{{\boldsymbol{p}}}^{\left({k}_{1},{k}_{2},\cdots ,{k}_{t}\right)},\end{eqnarray}$in which ${{T}}_{{k}_{i}}=\left(2{k}_{i}-1\right)!!{J}_{{E}_{{k}_{i}},{F}_{{k}_{i}}}$ and the operator ${{T}}_{{k}_{i}}$ just acts on those indices of the corresponding rank k_t of ${{\boldsymbol{p}}}^{\left({k}_{1},{k}_{2},\cdots ,{k}_{t}\right)}.$ When t=1, ${{\boldsymbol{\xi }}}^{\left({k}_{1}\right)}$ is the traceless multipole with ξ⁽¹⁾=μ, ξ⁽²⁾=Θ, ξ⁽³⁾=Ω, ξ⁽⁴⁾=Φ, ···. With the help of equations (C5), (C6), (C11) and (C12), the relationship between Applequist’s and Buckingham’s traceless multipoles or (hyper)polarizabilities is:(C13)$ \begin{eqnarray}\begin{array}{lll}{{\boldsymbol{\xi }}}^{\left({k}_{1},{k}_{2},\cdots ,{k}_{t}\right)} & = & \displaystyle \prod _{M}^{\left\{M| {n}_{M}\ne 0,M=2,3,\cdots \right\}}{n}_{M}!{\left({f}_{M}\right)}^{{n}_{M}-1}{{\boldsymbol{\Xi }}}_{{n}_{1}{n}_{2}{n}_{3}\cdots }\\ & = & \displaystyle \prod _{M}^{\left\{M| {n}_{M}\ne 0,M=2,3,\cdots \right\}}{n}_{M}!{\left[\displaystyle \frac{\left(2M-1\right)!!}{M!}\right]}^{{n}_{M}-1}{{\boldsymbol{\Xi }}}_{{n}_{1}{n}_{2}{n}_{3}\cdots },\end{array}\end{eqnarray}$in which ${{\rm{\Xi }}}_{{n}_{1}{n}_{2}{n}_{3}\cdots }$ is the traceless multipole or (hyper)polarizability defined in the equation (C5) and n₁+n₂+ n₃+ ···=t.

Reference View Option By original order
By published year
By cited within times
By Impact factor

[1]Buckingham A D Fowler P W Hutson J M 1988 Theoretical studies of van der Waals molecules and intermolecular forces
Chem. Rev. 88 963 988
DOI:10.1021/cr00088a008 [Cited within: 1]

[2]Stone A J 2013 The Theory of Intermolecular Forces Oxford Oxford University Press
[Cited within: 5]

[3]Buckingham A D 1959 Direct method of measuring molecular quadrupole moments
J. Chem. Phys. 30 1580 1585
DOI:10.1063/1.1730242 [Cited within: 9]

[4]Buckingham A D 1959 Molecular quadrupole moments
Quart. Rev. Chem. Soc. 13 183 214
DOI:10.1039/qr9591300183 [Cited within: 3]

[5]Buckingham A D 1967 Permanent and induced molecular moments and long-range intermolecular forces
Adv. Chem. Phys. 12 107 142
DOI:10.1002/9780470143582.ch2 [Cited within: 6]

[6]Buckingham A D Pople J A 1955 Theoretical studies of the kerr effect: I. Deviations from a linear polarization law
Proc. Phys. Soc. A 68 905 909
DOI:10.1088/0370-1298/68/10/307 [Cited within: 3]

[7]Buckingham A D Disch R L 1963 The quadrupole moment of the carbon dioxide molecule
Proc. R. Soc. A 273 275 289
DOI:10.1098/rspa.1963.0088 [Cited within: 4]

[8]Bishop D M Pipin J 1987 Field and field-gradient polarizabilities of H₂O
Theor. Chim. Acta 71 247 253
DOI:10.1007/BF00529096 [Cited within: 3]

[9]Maroulis G 1998 Hyperpolarizability of H₂O revisited: accurate estimate of the basis set limit and the size of electron correlation effects
Chem. Phys. Lett. 289 403 411
DOI:10.1016/S0009-2614(98)00439-4 [Cited within: 2]

[10]Haskopoulos A Maroulis G 2006 Dipole and quadrupole (hyper) polarizability for the asymmetric stretching of carbon dioxide: improved agreement between theory and experiment
Chem. Phys. Lett. 417 235 240
DOI:10.1016/j.cplett.2005.10.023 [Cited within: 4]

[11]Elking D M Perera L Duke R Darden T Pedersen L G 2011 A finite field method for calculating molecular polarizability tensors for arbitrary multipole rank
J. Comput. Chem. 32 3283 3295
DOI:10.1002/jcc.21914 [Cited within: 4]

[12]McLean A D Yoshimine M 1967 Theory of molecular polarizabilities
J. Chem. Phys. 47 1927 1935
DOI:10.1063/1.1712220 [Cited within: 18]

[13]Applequist J 1984 Fundamental relationships in the theory of electric multipole moments and multipole polarizabilities in static fields
Chem. Phys. 85 279 290
DOI:10.1016/0301-0104(84)85039-9 [Cited within: 6]

[14]Thakkar A J Lupinetti C 2006 Atomic polarizabilities and hyperpolarizabilities: a critical compilation
Atoms, Molecules and Clusters in Electric Fields: Theoretical Approaches to the Calculation of Electric Polarizability.Maroulis GLondon Imperial College Press
[Cited within: 1]

[15]Frisch M J et al. 2016 Gaussian 16, Revision B.01. Wallingford, CT Gaussian, Inc.
[Cited within: 3]

[16]Rhoderick E H 1947 Definition of nuclear quadrupole moments
Nature 160 255 256
DOI:10.1038/160255b0 [Cited within: 1]

[17]Chetty N Couling V W 2011 Measurement of the electric quadrupole moments of CO₂ and OCS
Mol. Phys. 109 655 666
DOI:10.1080/00268976.2010.546375 [Cited within: 2]

[18]Watson J N Craven I E Ritchie G L D 1997 Temperature dependence of electric field-gradient induced birefringence in carbon dioxide and carbon disulfide
Chem. Phys. Lett. 274 1 6
DOI:10.1016/S0009-2614(97)00656-8 [Cited within: 1]

[19]Graham C Imrie D A Raab R E 1998 Measurement of the electric quadrupole moments of CO₂, CO, N₂, Cl₂ and BF₃
Mol. Phys. 93 49 56
DOI:10.1080/00268979809482187 [Cited within: 2]

[20]Balachandran Pillai P C Couling V W 2019 Dispersion of the Rayleigh light-scattering virial coefficients and polarisability anisotropy of CO₂
Mol. Phys. 117 289 297
DOI:10.1080/00268976.2018.1510143 [Cited within: 2]

[21]Gentle I R Laver D R Ritchie G L D 1989 Second hyperpolarizability and static polarizability anisotropy of carbon dioxide
J. Phys. Chem. 93 3035 3038
DOI:10.1021/j100345a033 [Cited within: 1]

[22]Maroulis G Thakkar A J 1990 Polarizabilities and hyperpolarizabilities of carbon dioxide
J. Chem. Phys. 93 4164 4171
DOI:10.1063/1.458749 [Cited within: 5]

[23]Maroulis G 2003 Electric (hyper) polarizability derivatives for the symmetric stretching of carbon dioxide
Chem. Phys. 291 81 95
DOI:10.1016/S0301-0104(03)00186-1 [Cited within: 5]

[24]Graner G Rossetti C Bailly D 1986 The carbon dioxide molecule: a test case for the r₀, r_e and r_m structures
Mol. Phys. 58 627 636
DOI:10.1080/00268978600101431 [Cited within: 1]

[25]Kurtz H A Stewart J J Dieter K M 1990 Calculation of the nonlinear optical properties of molecules
J. Comput. Chem. 11 82 87
DOI:10.1002/jcc.540110110 [Cited within: 2]

[26]Alms G R Burnham A K Flygare W H 1975 Measurement of the dispersion in polarizability anisotropies
J. Chem. Phys. 63 3321 3326
DOI:10.1063/1.431821 [Cited within: 1]

[27]Bogaard M P Buckingham A D Pierens R K White A H 1978 Rayleigh scattering depolarization ratio and molecular polarizability anisotropy for gases
J. Chem. Soc. Faraday Trans. 1 74 3008 3015
DOI:10.1039/f19787403008 [Cited within: 1]

[28]Maroulis G 2004 A note on the electric multipole moments of carbon dioxide
Chem. Phys. Lett. 396 66 68
DOI:10.1016/j.cplett.2004.07.123 [Cited within: 1]

[29]Loboda O Ingrosso F Ruiz-López M F Reis H Millot C 2016 Dipole and quadrupole polarizabilities of the water molecule as a function of geometry
J. Comput. Chem. 37 2125 2132
DOI:10.1002/jcc.24431 [Cited within: 1]

[30]Nénon S Champagne B Spassova M I 2014 Assessing long-range corrected functionals with physically-adjusted range-separated parameters for calculating the polarizability and the second hyperpolarizability of polydiacetylene and polybutatriene chains
Phys. Chem. Chem. Phys. 16 7083 7088
DOI:10.1039/c4cp00105b [Cited within: 1]

[31]Oviedo M B Ilawe N V Wong B M 2016 Polarizabilities of π-conjugated chains revisited: improved results from broken-symmetry range-separated DFT and new CCSD(T) benchmarks
J. Chem. Theory Comput. 12 3593 3602
DOI:10.1021/acs.jctc.6b00360 [Cited within: 1]

[32]Burgos E Bonadeo H 1981 Electrical multipoles and multipole interactions: compact expressions and a diagrammatic method
Mol. Phys. 44 1 15
DOI:10.1080/00268978100102251 [Cited within: 1]