ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/multipole/dielectric_new.tex
Revision: 4353
Committed: Tue Sep 15 15:56:07 2015 UTC (10 years ago) by mlamichh
Content type: application/x-tex
File size: 48885 byte(s)
Log Message:
Added new tex file for dielectric.

File Contents

# User Rev Content
1 mlamichh 4353 \documentclass[%
2     aip,
3     jcp,
4     amsmath,amssymb,
5     preprint,%
6     % reprint,%
7     %author-year,%
8     %author-numerical,%
9     ]{revtex4-1}
10    
11     \usepackage{graphicx}% Include figure files
12     \usepackage{dcolumn}% Align table columns on decimal point
13     \usepackage{multirow}
14     \usepackage{bm}% bold math
15     \usepackage{natbib}
16     \usepackage{times}
17     \usepackage{mathptmx}
18     \usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions
19     \usepackage{url}
20     \usepackage{braket}
21    
22     \begin{document}
23    
24     \title{Real space electrostatics for multipoles. III. Dielectric Properties}
25    
26     \author{Madan Lamichhane}
27     \affiliation{Department of Physics, University
28     of Notre Dame, Notre Dame, IN 46556}
29     \author{Thomas Parsons}
30     \affiliation{Department of Chemistry and Biochemistry, University
31     of Notre Dame, Notre Dame, IN 46556}
32     \author{Kathie E. Newman}
33     \affiliation{Department of Physics, University
34     of Notre Dame, Notre Dame, IN 46556}
35     \author{J. Daniel Gezelter}
36     \email{gezelter@nd.edu.}
37     \affiliation{Department of Chemistry and Biochemistry, University
38     of Notre Dame, Notre Dame, IN 46556}
39    
40     \date{\today}% It is always \today, today,
41     % but any date may be explicitly specified
42    
43     \begin{abstract}
44     Note: This manuscript is a work in progress.
45    
46     We report on the dielectric properties of the shifted potential
47     (SP), gradient shifted force (GSF), and Taylor shifted force (TSF)
48     real-space methods for multipole interactions that were developed in
49     the first two papers in this series. We find that some subtlety is
50     required for computing dielectric properties with the real-space
51     methods, particularly when using the common fluctuation formulae.
52     Three distinct methods for computing the dielectric constant are
53     investigated, including the standard fluctuation formulae,
54     potentials of mean force between solvated ions, and direct
55     measurement of linear solvent polarization in response to applied
56     fields and field gradients.
57     \end{abstract}
58    
59     \maketitle
60    
61     \section{Introduction}
62    
63     Over the past several years, there has been increasing interest in
64     pairwise methods for correcting electrostatic interactions in computer
65     simulations of condensed molecular
66     systems.\cite{Wolf99,Zahn02,Kast03,Beckd.A.C._Bi0486381,Ma05,Fennell06}
67     These techniques were initially developed from the observations and efforts of
68     Wolf {\it et al.} and their work towards an $\mathcal{O}(N)$
69     Coulombic sum.\cite{Wolf99} Wolf's method of cutoff neutralization is
70     able to obtain excellent agreement with Madelung energies in ionic
71     crystals.\cite{Wolf99} Later, Zahn \textit{et al.} and Fennell and Gezelter extended this method which incorporates Wolf's electrostatic energy and modified it to conserve the total energy in molecular dynamic simulation.\cite{Zahn02, Fennell06} We have also previously developed three generalized real space methods: Shifted potential (SP), Gradeint shifted force (GSF), and Taylor shifted force (TSF) for evaluating electrostatic interactions for higher order multipoles (dipoles and quadrupoles) in the previous two papers in the series.\cite{PaperI, PaperII} These real space methods consider a finite cutoff sphere with the neutralization of the electrostatic moment within the cutoff sphere. Furthermore, extra terms have been added to the potential energy to vanish force and torque smoothly at the cutoff radius which ensures the conservation of the total energy in the molecular dynamic simulation.
72    
73     One of the most difficult tests of any new electrostatic method is the fidelity with which that method can reproduce the bulk-phase polarizability or equivalently, the dielectric properties of a fluid. Since dielectric properties are macroscopic properties, all interactions between molecules in an entire system are significantly important. But it is computationally infeasible to consider every interactions between molecules in the macroscopically large system. Therefore small molecular system with periodic boundary condition and finite cutoff region of interactions is usually considered in computer simulations. While calculating dielectric properties, the formula should be modified in such a way so that it can accommodate behaviour of electrostatic neutrality and smoothness of energy, force and torque at the cutoff radius. Previously many studies have been conducted to calculate dipolar and quadrupolar dielectric properties using computer simulations. \cite{Kirkwood39, Onsagar36,LoganI81, LoganII82, LoganIII82} But these methods do not specifically take account of the cutoff behavior common in real-space electrosatic methods. In 1983 Neumann proposed a general formula for evaluating dielectric properties for dipolar fluid using real-space cutoff methods. \cite{Neumann83} In the same year Steinhauser and Neumann used this formula to evaluate the correct dielectric constant for the Stockmayer fluid using two different methods: Ewald-Kornfield (EK) and reaction field (RF) methods. \cite{Neumann-Steinhauser83} This formula contains a correction factor which is equal to $\frac{3}{4 \pi} $ times volume integral of the dipole-dipole interactions for a given electrostatic cutoff method (See equation \ref{dipole-diopleTensor}).\cite{Neumann83} Similarly Zahn \textit{et al.}\cite{Zahn02} also evaluated correction factor for dipole-dipole interaction using damped shifted charge-charge kernel (see equation \ref{dipole-chargeTensor}). This later generalized by Izvekove \textit{et al.}, which is equal to $\frac{3}{4 \pi} $. \cite{Izvekov:2008wo} When the correction factor is equal to $\frac{3}{4 \pi} $, the expression for the dielectric constant reduces to widely-used \textit{conducting boundary} formula (see equation (\ref{correctionFormula})). Many studies have also been conducted to understand solvation theory using dielectric properties of quadrupolar fluid.\cite{JeonI03, JeonII03, Chitanvis96}. But these studies do not use correction factor straight forwardly to evaluate correct dielectric properties for quadrupolar fluid.
74    
75     In this paper we are proposing analogous general formula for calculating dielectric properties for quadrupolar fluid. Furthermore we have also evaluated the correction factor for SP, GSF, and TSF method for both dipolar and quadrupolar fluid considering charge-charge, dipole-dipole or quadrupole-quadrupole interactions. The relation between quadrupolar susceptibility and dielectric constant is not straight forward for quadrupolar fluid as in the dipolar case. The dielectric constant depends on the geometry of the external field perturbation.\cite{Ernst92} We have also calculated the geometrical factor for two ions immersed quadrupolar system to evaluate dielectric constant from the quadrupolar susceptibility. We have used three different methods: fluctuation formula, external field perturbation and the potential of mean force, to study dielectric properties of the dipolar and quadrupolar system. The fluctuation formula observes the time average fluctuation of the multipolar moment as a function of temperature. The average fluctuation value of the system is determined by the multipole-multipole interactions between molecules at a given temperature. In the external field perturbation, the net polarization of the system is observed as a linear response of the applied field perturbation, where proportionality constant is determined by the electrostatic interaction between the electrostatic multipoles at a given temperature. Since the expression of the electrostatic interaction energy, force, and torque in the real space electrostatic methods are slightly modified from their original definition therefore both fluctuation and external field perturbation formula should also be modified accordingly.The potential of mean force method calculates dielectric constant from the potential energy between ions before and after dielectric material is introduced. All of these different methods for calculating dielectric properties will be discussed in detail in the following sections: \ref{sec:perturbation}, \ref{sec:fluctuation}, and \ref{sec:PMF}.
76     \subsection{Boltzmann average of orientational polarization - dipole}
77     Consider a system of molecules with permenent dipole moment $p_o$. In the absense of external field, thermal agitation makes dipole randomly oriented therefore there is no net dipole moment. But external field tends them to line up in the direction of applied field. Therefore the total Hamiltonian of the molecule is,\cite{Jackson98}
78    
79     \begin{equation}
80     H = H_o - \bf{p_o} .\bf{E},
81     \end{equation}
82     where $H_o$ is a function of the internal coordinates of the molecule. Now Boltzmann average of the dipole moment is given by,
83     \begin{equation}
84     \braket{p_{mol}} = \frac{\displaystyle\int d\Omega\; p_o\; cos\theta\; e^{\frac{p_oE\; cos\theta}{k_B T}}}{\displaystyle\int d\Omega\; e^{\frac{p_oE\;cos\theta}{k_B T}}},
85     \end{equation}
86     where $\bf{E}$ is selected along z-axis. If we consider applied field is small i.e. $\frac{p_oE\; cos\theta}{k_B T} << 1$ then we get,
87    
88     \begin{equation}
89     \braket{p_{mol}} \approx \frac{1}{3}\frac{{p_o}^2}{k_B T}E,
90     \end{equation}
91     where $ \alpha_p = \frac{1}{3}\frac{{p_o}^2}{k_B T}$ is a molecular polarizability. The orientational polarization depends inversely on the temperature and applied field must overcome the thermal agitation.
92     \label{boltzAverage-dipole}
93    
94     \subsection{Boltzmann average of orientational polarization - quadrupole}
95     Consider a system of molecules with permanent quadrupole moment $q_{\alpha\beta} $. The average quadrupole moment at temperature T in the presence of uniform applied field gradient is given by,\cite{AduGyamfi78, AduGyamfi81}
96     \begin{equation}
97     \braket{q_{\alpha\beta}} \;=\; \frac{\displaystyle\int d\Omega\; e^{-\frac{H}{k_B T}}q_{\alpha\beta}}{\displaystyle\int d\Omega\; e^{-\frac{H}{k_B T}}} \;=\; \frac{\displaystyle\int d\Omega\; e^{\frac{q_{\mu\nu}\;\partial_\nu E_\mu}{k_B T}}q_{\alpha\beta}}{\displaystyle\int d\Omega\; e^{\frac{q_{\mu\nu}\;\partial_\nu E_\mu}{k_B T}}},
98     \label{boltzQuad}
99     \end{equation}
100     where $\int d\Omega = \int_0^{2\pi} \int_0^\pi \int_0^{2\pi}
101     sin\theta\; d\theta\ d\phi\ d\psi$ is the integration over Euler
102     angles, $ H = H_o -q_{\mu\nu}\;\partial_\nu E_\mu $ is the energy of
103     a quadrupole in the gradient of the
104     applied field and $ H_o$ is a function of internal coordinates of the molecule. The energy and quadrupole moment can be transformed into body frame using following relation,
105     \begin{equation}
106     \begin{split}
107     &q_{\alpha\beta} = \eta_{\alpha\alpha'}{q}^* _{\alpha'\beta'}\eta_{\beta'\beta} \\
108     &H = H_o - q:\vec{\nabla}\vec{E} = H_o - q_{\mu\nu}\;\partial_\nu E_\mu = H_o -\eta_{\mu\mu'}{q}^*_{\mu'\nu'}\eta_{\nu'\nu}\;\partial_\nu E_\mu.
109     \end{split}
110     \label{energyQuad}
111     \end{equation}
112     Here the starred tensors are the components in the body fixed
113     frame. Substituting equation (\ref{energyQuad}) in the equation (\ref{boltzQuad})
114     and taking linear terms in the expansion we get,
115     \begin{equation}
116     \braket{q_{\alpha\beta}} = \frac{ \int d\Omega \left(1 + \frac{\eta_{\mu\mu'}{q}^*_{\mu'\nu'}\eta_{\nu'\nu}\;\partial_\nu E_\mu }{k_B T}\right)q_{\alpha\beta}}{ \int d\Omega \left(1 + \frac{\eta_{\mu\mu'}{q}^*_{\mu'\nu'}\eta_{\nu'\nu}\;\partial_\nu E_\mu }{k_B T}\right)},
117     \end{equation}
118     where $\eta_{\alpha\alpha'}$ is the rotation matrix that transforms
119     the body fixed to the space fixed frame,
120     \[\eta_{\alpha\alpha'}
121     = \left(\begin{array}{ccc}
122     cos\phi\; cos\psi - cos\theta\; sin\phi\; sin\psi & cos\theta\; cos\psi\; sin\phi + cos\phi\; sin\psi & sin\theta\; sin\phi \\
123     -cos\psi\; sin\phi - cos\theta\; cos\phi \; sin\psi & cos\theta\; cos\phi\; cos\psi - sin\phi\; sin\psi & cos\phi\; cos\theta \\
124     sin\theta\; sin\psi & -cos\psi\; sin\theta & cos\theta
125     \end{array} \right).\]
126     Integration of 1st and 2nd terms in the denominator gives $8 \pi^2$
127     and $8 \pi^2 /3\;\vec{\nabla}.\vec{E}\; Tr(q^*) $ respectively. The
128     second term vanishes for charge free space
129     (i.e. $\vec{\nabla}.\vec{E} \; = \; 0)$. Similarly integration of the
130     1st term in the numerator produces
131     $8 \pi^2 /3\; Tr(q^*)\delta_{\alpha\beta}$ and the 2nd term produces
132     $8 \pi^2 /15k_B T (3{q}^*_{\alpha'\beta'}{q}^*_{\beta'\alpha'} -
133     {q}^*_{\alpha'\alpha'}{q}^*_{\beta'\beta'})\partial_\alpha E_\beta$,
134     if $\vec{\nabla}.\vec{E} \; = \; 0$,
135     $ \partial_\alpha E_\beta = \partial_\beta E_\alpha$ and
136     ${q}^*_{\alpha'\beta'}= {q}^*_{\beta'\alpha'}$. Therefore the
137     Boltzmann average of a quadrupole moment can be written as,
138    
139     \begin{equation}
140     \braket{q_{\alpha\beta}}\; = \; \frac{1}{3} Tr(q^*)\;\delta_{\alpha\beta} + \frac{{\bar{q_o}}^2}{15k_BT}\;\partial_\alpha E_\beta,
141     \end{equation}
142     where $ \alpha_q = \frac{{\bar{q_o}}^2}{15k_BT} $ is a molecular quadrupolarizablity and ${\bar{q_o}}^2=
143     3{q}^*_{\alpha'\beta'}{q}^*_{\beta'\alpha'}-{q}^*_{\alpha'\alpha'}{q}^*_{\beta'\beta'}$ is a square of the net quadrupole moment of a molecule.
144    
145     \section{External field perturbation}
146     In the presence of uniform electric field $\textbf{E}^o$ (or field gradient $\partial_\alpha E^o_\beta$), a system of dipolar (or quadrupolar) molecules polarizes along the direction of the applied field (or field gradient). Therefore the net dipolar polarization $ \textbf{P}$ of the system is,
147     \begin{equation}
148     \textbf{P} = \epsilon_o \alpha_{D}\; \textbf{E}^o.
149     \label{pertDipole}
150     \end{equation}
151     The constant $\alpha_D$ is a macroscopic polarizability, which is a property of the dipolar fluid in a given density and temperature. Similarly, the net quadrupolar polarization of the system can be given by,
152     \begin{equation}
153     \begin{split}
154     & {Q}_{\alpha\beta} = \frac{1}{3}\; Tr({Q})\; \delta_{\alpha\beta} + \epsilon_o\; \alpha_Q \; \partial_{\alpha} E^o_{\beta}
155     \\ & or \\
156     & \frac{1}{3}\;\Theta_{\alpha\beta} = \epsilon_o\; \alpha_Q \; \partial_{\alpha} E^o_{\beta}
157     \end{split}
158     \label{pertQuad}
159     \end{equation}
160     where $Q_{\alpha\beta}$ is a tensor component of the traced quadrupolar moment of the system, $ \alpha_Q$ is a macroscopic quadrupolarizability has a dimension of $length^{-2}$, and $\Theta_{\alpha\beta} = 3Q_{\alpha\beta}-Tr(Q) $ is the traceless component of the quadrupole moment.
161     \label{sec:perturbation}
162    
163     \section{Fluctuation formula}
164     For a system of molecules with net dipolar moment $\bf{M}$ at thermal equilibrium of temperature T in the presence of applied field $\bf{E}^o$, the average dipolar polarization can be expressed in terms of fluctuation of the net dipole moment as below,\cite{Stern03}
165     \begin{equation}
166     \braket{\bf{P}} = \epsilon_o \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T}\bf{E}^o
167     \label{flucDipole}
168     \end{equation}
169     This is similar to the formula for boltzmann average of single dipolar molecule in the subsection \ref{boltzAverage-dipole}. Here $\braket{\bf{P}}$ is average polarization and $ \braket{\textbf{M}^2}-{\braket{\textbf{M}}}^2$ is the net dipole fluctuation at temperature T. For the limiting case $\textbf{E}^o \rightarrow 0 $, ensemble average of both net dipole moment $\braket{\textbf{M}}$ and dipolar polarization $\braket{\bf{P}}$ tends to vanish but $\braket{\bf{M}^2}$ will still be non-zero. The dipolar macroscopic polarizability can be written as,
170     \begin{equation}
171     \alpha_D = \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T}
172     \end{equation}
173     This is a macroscopic property of dipolar material which is true even if applied field $ \textbf{E}^o \rightarrow 0 $.
174    
175     Analogous formula can also be written for a system with quadrupolar molecules,
176     \begin{equation}
177     \braket{Q_{\alpha\beta}} = \frac{1}{3} Tr(\textbf{Q})\; \delta_{\alpha\beta} + \epsilon_o \frac{\braket{\textbf{Q}^2}-{\braket{\textbf{Q}}}^2}{15 \epsilon_o V k_B T}{\partial_\alpha E^o_\beta}
178     \label{flucQuad}
179     \end{equation}
180     where $Q_{\alpha\beta}$ is a component of system quadrupole moment, $\bf{Q}$ is net quadrupolar moment which can be expressed as $\textbf{Q}^2 =3Q_{\alpha\beta}Q_{\alpha\beta}-(Tr\textbf{Q})^2 $. The macroscopic quadrupolarizability is given by,
181     \begin{equation}
182     \alpha_Q = \frac{\braket{\textbf{Q}^2}-{\braket{\textbf{Q}}}^2}{15 \epsilon_o V k_B T}
183     \label{propConstQuad}
184     \end{equation}
185     \label{sec:fluctuation}
186    
187     \section{Potential of mean force}
188     In this method, we will measure the interaction between a positive and negative charge at varying distances after introducing a dipolar (or quadrupolar) material between them. The potential of mean force (PMF) between two ions in a liquid is obtained by constraining their distance and measuring the mean constraint force required to hold them at a fixed distance $r.$ The PMF is obtained from a sequence of simulations as,
189     \begin{equation}
190     w(r) = \int_{\inf}^{r}\braket{\frac{\partial f}{\partial r'}}dr',
191     \end{equation}
192     where $\braket{\partial f/\partial r'}$ is the mean constraint force.
193     Since the ions have a protecting Lennard-Jones (LJ) potential,
194     \begin{equation}
195     w(r) = w_{LJ}(r) + \frac{q_iq_j}{4\pi \epsilon_o \epsilon(r)}U_{method}(r).
196     \end{equation}
197     Here $w_{LJ}$ is the PMF calculated without electrostatic interactions and $U_{method}(r)$ is the radial function for the charge-charge interaction, which is different for various real space truncation methods.
198    
199     The quadrupole molecule can only couple with the gradient of the electric field and the region between two opposite point charges has both an electric field and a gradient of the electric field present. Therefore, this methodology should be usable to determine the dielectric constant for both the dipolar and quadrupolar fluid.
200     \label{sec:PMF}
201    
202     \section{Correction factor}
203     Since equations (\ref{pertDipole}, \ref{pertQuad}, \ref{flucDipole}, and \ref{flucQuad}) provide relation between polarization (dipolar or quadrupolar) and applied field (uniform field or field gradient), $\chi_d$ (or $ \chi_q$) is actually a macroscopic polarizability (or quadrupolarizability), which is different than the dipolar (or quadrupolar) susceptibility of the fluid. Actual constitutive relation should have a relation between polarization and Maxwell field (or field gradient) at different point in the sample. We can obtain susceptibility of the fluid from its macroscopic polarizability using correction factor evaluated below.
204     \subsection{Dipolar system}
205     In the presence of an external field $ \textbf{E}$ polarization $\textbf{E}$ will be induced in a dipolar system. The total electrostatic field (or Maxwell electric field) at point $\bf{r}$ in a system is,\cite{Neumann83}
206     \begin{equation}
207     \textbf{E}(\textbf{r}) = \textbf{E}^o(\textbf{r}) + \frac{1}{4\pi\epsilon_o} \int d^3r' \textbf{T}(\textbf{r}-\textbf{r}')\cdot {\textbf{P}(\textbf{r}')}.
208     \end{equation}
209    
210     We can consider the cases of Stockmayer (dipolar) soft spheres that are represented either by two closely-spaced point charges or by a single point dipole (see Fig. \ref{fig:stockmayer}).
211     \begin{figure}
212     \includegraphics[width=3in]{DielectricFigure}
213     \caption{With the real-space electrostatic methods, the effective
214     dipole tensor, $\mathbf{T}$, governing interactions between
215     molecular dipoles is not the same for charge-charge interactions as
216     for point dipoles.}
217     \label{fig:stockmayer}
218     \end{figure}
219     In the case where point charges are interacting via an electrostatic
220     kernel, $v(r)$, the effective {\it molecular} dipole tensor,
221     $\mathbf{T}$ is obtained from two successive applications of the
222     gradient operator to the electrostatic kernel,
223     \begin{equation}
224     \mathbf{T}_{\alpha \beta}(r) = \nabla_\alpha \nabla_\beta \left(v(r)\right) = \delta_{\alpha \beta}
225     \left(\frac{1}{r} v^\prime(r) \right) + \frac{r_{\alpha}
226     r_{\beta}}{r^2} \left( v^{\prime \prime}(r) - \frac{1}{r}
227     v^{\prime}(r) \right)
228     \label{dipole-chargeTensor}
229     \end{equation}
230     where $v(r)$ may be either the bare kernel ($1/r$) or one of the
231     modified (Wolf or DSF) kernels. This tensor describes the effective
232     interaction between molecular dipoles ($\mathbf{D}$) in Gaussian
233     units as $-\mathbf{D} \cdot \mathbf{T} \cdot \mathbf{D}$.
234    
235     When utilizing the new real-space methods for point dipoles, the
236     tensor is explicitly constructed,
237     \begin{equation}
238     \mathbf{T}_{\alpha \beta}(r) = \delta_{\alpha \beta} v_{21}(r) +
239     \frac{r_{\alpha} r_{\beta}}{r^2} v_{22}(r)
240     \label{dipole-diopleTensor}
241     \end{equation}
242     where the functions $v_{21}(r)$ and $v_{22}(r)$ depend on the level of
243     the approximation. Although the Taylor-shifted (TSF) and
244     gradient-shifted (GSF) models produce to the same $v(r)$ function for
245     point charges, they have distinct forms for the dipole-dipole
246     interactions.
247    
248     Using constitutive relation, the polarization density $\textbf{P}(\textbf{r})$ is given by,
249     \begin{equation}
250     \textbf{P}(\textbf{r}) = \epsilon_o\; \chi^*_D \left(\textbf{E}^o(\textbf{r}) + \frac{1}{4\pi\epsilon_o} \int d^3r' \textbf{T}(\textbf{r}-\textbf{r}')\cdot {\textbf{P}(\textbf{r}')}\right).
251     \label{constDipole}
252     \end{equation}
253     Here $\chi^*_D$ is a dipolar susceptibility can be expressed in terms of dielectric constant as $ \chi^*_D = \epsilon - 1$ which different than macroscopic dipolar polarizability $\alpha_D$ in the sections \ref{sec:perturbation} and \ref{sec:fluctuation}. We can split integral into two parts: singular part i.e $|\textbf{r}-\textbf{r}'|\rightarrow 0 $ and non-singular part i.e $|\textbf{r}-\textbf{r}'| > 0 $ . The singular part of the integral can be written as,\cite{Neumann83, Jackson98}
254     \begin{equation}
255     \frac{1}{4\pi\epsilon_o} \int_{|\textbf{r}-\textbf{r}'| \rightarrow 0} d^3r'\; \textbf{T}(\textbf{r}-\textbf{r}')\cdot {\textbf{P}(\textbf{r}')} = - \frac{\textbf{P}(\textbf{r})}{3\epsilon_o}
256     \label{singular}
257     \end{equation}
258     Substituting equation (\ref{singular}) in the equation (\ref{constDipole}) we get,
259     \begin{equation}
260     \textbf{P}(\textbf{r}) = 3 \epsilon_o\; \frac{\chi^*_D}{\chi^*_D + 3} \left(\textbf{E}^o(\textbf{r}) + \frac{1}{4\pi\epsilon_o} \int_{|\textbf{r}-\textbf{r}'| > 0} d^3r'\; \textbf{T}(\textbf{r}-\textbf{r}')\cdot {\textbf{P}(\textbf{r}')}\right).
261     \end{equation}
262     For both polarization and electric field homogeneous, this can be easily solved using Fourier transformation,
263     \begin{equation}
264     \textbf{P}(\kappa) = 3 \epsilon_o\; \frac{\chi^*_D}{\chi^*_D + 3} \left(1- \frac{3}{4\pi}\;\frac{\chi^*_D}{\chi^*_D + 3}\; \textbf{T}({\kappa})\right)^{-1}\textbf{E}^o({\kappa}).
265     \end{equation}
266     For homogeneous applied field Fourier component is non-zero only if $\kappa = 0$. Therefore,
267     \begin{equation}
268     \textbf{P}(0) = 3 \epsilon_o\; \frac{\chi^*_D}{\chi^*_D + 3} \left(1- \frac{\chi^*_D}{\chi^*_D + 3}\; A_{dipole})\right)^{-1}\textbf{E}^o({0}).
269     \label{fourierDipole}
270     \end{equation}
271     where $A_{dipole}=\frac{3}{4\pi}T(0) = \frac{3}{4\pi} \int_V d^3r\;T(r)$. Now equation (\ref{fourierDipole}) can be compared with equation (\ref{flucDipole}). Therefore,
272     \begin{equation}
273     \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T} = \frac{3\;\chi^*_D}{\chi^*_D + 3} \left(1- \frac{\chi^*_D}{\chi^*_D + 3}\; A_{dipole})\right)^{-1}
274     \end{equation}
275     Substituting $\chi^*_D = \epsilon-1$ and $ \frac{\braket{\bf{M}^2}-{\braket{\bf{M}}}^2}{3 \epsilon_o V k_B T} = \epsilon_{CB}-1$ in above equation we get,
276     \begin{equation}
277     \epsilon = \frac{3+(A_{dipole} + 2)(\epsilon_{CB}-1)}{3+(A_{dipole} -1)(\epsilon_{CB}-1)}
278     \label{correctionFormula}
279     \end{equation}
280     where $\epsilon_{CB}$ is dielectric constant obtained from conducting boundary condition. Equation (\ref{correctionFormula}) calculates actual dielectric constant from the dielectric constant obtained from the conducting boundary condition (which can be obtained directly from the simulation) using correction factor. The correction factor is different for different real-space cutoff methods. The expression for correction factor assuming a single point dipole or two closely spaced point charges for SP, GSF, and TSF method is listed in Table \ref{tab:A}.
281     \begin{table}
282     \caption{Expressions for the dipolar correction factor ($A$) for the real-space electrostatic methods in terms of the damping parameter
283     ($\alpha$) and the cutoff radius ($r_c$). The Ewald-Kornfeld result
284     derived in Refs. \onlinecite{Adams:1980rt,Adams:1981fr,Neumann83} is shown for comparison using the Ewald
285     convergence parameter ($\kappa$) and the real-space cutoff value ($r_c$). }
286     \label{tab:A}
287     {%
288     \begin{tabular}{l|c|c|c|}
289    
290     Method & $A_\mathrm{charges}$ & $A_\mathrm{dipoles}$ \\
291     \hline
292     Spherical Cutoff (SC) & $\mathrm{erf}(r_c \alpha) - \frac{2 \alpha r_c}{\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $\mathrm{erf}(r_c \alpha) - \frac{2 \alpha r_c}{\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ \\
293     Shifted Potental (SP) & $ \mathrm{erf}(r_c \alpha) - \frac{2 \alpha r_c}{\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $\mathrm{erf}(r_c \alpha) -\frac{2 \alpha r_c}{\sqrt{\pi}}\left(1+\frac{2\alpha^2 {r_c}^2}{3} \right)e^{-\alpha^2{r_c}^2} $\\
294     Gradient-shifted (GSF) & 1 & $\mathrm{erf}(\alpha r_c)-\frac{2 \alpha r_c}{\sqrt{\pi}} \left(1 + \frac{2 \alpha^2 r_c^2}{3} + \frac{\alpha^4 r_c^4}{3}\right)e^{-\alpha^2 r_c^2} $ \\
295     Taylor-shifted (TSF) & 1 & 1 \\
296     Ewald-Kornfeld (EK) & $\mathrm{erf}(r_c \kappa) - \frac{2 \kappa r_c}{\sqrt{\pi}} e^{-\kappa^2 r_c^2}$ & $\mathrm{erf}(r_c \kappa) - \frac{2 \kappa r_c}{\sqrt{\pi}} e^{-\kappa^2 r_c^2}$ \\\hline
297     \end{tabular}%
298     }
299     \end{table}
300     \subsection{Quadrupolar system}
301     In the presence of the field gradient $\partial_\alpha {E}_\beta $, a
302     non-vanishing quadrupolar polarization (quadrupole moment per unit
303     volume) $\bar{Q}_{\alpha\beta}$ will be induced in the entire volume
304     of a sample. The total field at any point $\vec{r}$ in the sample is
305     given by,
306     \begin{equation}
307     \partial_\alpha E_\beta(\textbf{r}) = \partial_\alpha {E^o}_\beta(\textbf{r}) + \frac{1}{4\pi \epsilon_o}\int T_{\alpha\beta\gamma\delta}(|{\textbf{r}-\textbf{r}'}|)\;{Q}_{\gamma\delta}(\textbf{r}')\; d^3r'
308     \label{gradMaxwell}
309     \end{equation}
310     where $\partial_\alpha {E^o}_\beta$ is the applied field gradient and $ T_{\alpha\beta\gamma\delta}$ is the quadrupole-quadrupole interaction tensor. We can represent quadrupole as a group of four closely spaced charges, two closely spaced point dipoles or single point quadrupole (see Fig. \ref{fig:quadrupolarFluid}). The quadrupole-quadrupole interaction tensor from the charge representation can obtained from the application of the four successive gradient operator to the electrostatic kernel $v(r)$.
311    
312     \begin{equation}
313     \begin{split}
314     T_{\alpha\beta\gamma\delta}(r) &=\nabla_\alpha \nabla_\beta \nabla_\gamma \nabla_\delta\;v(r)
315     \\ &= \left(\delta_{\alpha\beta}\delta_{\gamma\delta} + \delta_{\alpha\gamma}\delta_{\beta\delta}+ \delta_{\alpha\delta}\delta_{\beta\gamma}\right)\left(-\frac{v'(r)}{r^3} + \frac{v''(r)}{r^2}\right)
316     \\ &+ \left(\delta_{\alpha\beta} r_\gamma r_\delta + 5 \; permutations \right) \left(\frac{3v'(r)}{r^5}-\frac{3v''(r)}{r^4} + \frac{v'''(r)}{r^3}\right)
317     \\ &+ r_\alpha r_\beta r_\gamma r_\delta\; \left(-\frac{15v'(r)}{r^7}+\frac{15v''(r)}{r^6}-\frac{6v'''(r)}{r^5} + \frac{v''''(r)}{r^4}\right),
318     \end{split}
319     \label{quadCharge}
320     \end{equation}
321     where $v(r)$ can either be electrostatic kernel for spherical truncation or one of the modified (Wolf or DSF) method. Similarly in point dipole representation the qaudrupole-quadrupole interaction tensor can be obtained from the applications of the two successive gradient in the dipole-dipole interaction tensor,
322    
323     \begin{equation}
324     \begin{split}
325     T_{\alpha\beta\gamma\delta}(r) &=\nabla_\alpha \nabla_\beta \;v_{\gamma\delta}(r)
326     \\ &= \delta_{\alpha\beta}\delta_{\gamma\delta} \frac{v'_{21}(r)}{r} + \left(\delta_{\alpha\gamma}\delta_{\beta\delta}+ \delta_{\alpha\delta}\delta_{\beta\gamma}\right)\frac{v_{22}(r)}{r^2}
327     \\ &+ \delta_{\gamma\delta} r_\alpha r_\beta \left(\frac{v''_{21}(r)}{r^2}-\frac{v'_{21}(r)}{r^3} \right)
328     \\ &+\left(\delta_{\alpha\beta} r_\gamma r_\delta + \delta_{\alpha\gamma} r_\beta r_\delta +\delta_{\alpha\delta} r_\gamma r_\beta + \delta_{\beta\gamma} r_\alpha r_\delta +\delta_{\beta\delta} r_\alpha r_\gamma \right) \left(\frac{v'_{22}(r)}{r^3}-\frac{2v_{22}(r)}{r^4}\right)
329     \\ &+ r_\alpha r_\beta r_\gamma r_\delta\; \left(\frac{v''_{22}(r)}{r^4}-\frac{5v'_{22}(r)}{r^5}+\frac{8v_{22}(r)}{r^6}\right),
330     \end{split}
331     \label{quadDip}
332     \end{equation}
333     where $v_{\gamma\delta}(r)$ is the electrostatic dipole-dipole interaction tensor, which is different for different electrostatic cut off methods. Similarly $v_{21}(r) \;and\; v_{22}(r)$ are the radial function for different real space cutoff methods defined in Paper I of the series.\cite{PaperI} Using point quadrupole representation the quadrupole-quadrupole interaction can be constructed as,
334     \begin{equation}
335     \begin{split}
336     T_{\alpha\beta\gamma\delta}(r) &= \left(\delta_{\alpha\beta}\delta_{\gamma\delta} + \delta_{\alpha\gamma}\delta_{\beta\delta}+ \delta_{\alpha\delta}\delta_{\beta\gamma}\right)v_{41}(r) + \delta_{\gamma\delta} r_\alpha r_\beta \frac{v_{42}(r)}{r^2} \\ &+ r_\alpha r_\beta r_\gamma r_\delta\; \left(\frac{v_{43}(r)}{r^4}\right),
337     \end{split}
338     \label{quadRadial}
339     \end{equation}
340     where $v_{41}(r),\; v_{42}(r), \; \text{and} \; v_{43}(r)$ are defined in Paper I of the series. \cite{PaperI} They have different functional forms for different electrostatic cutoff methods.
341     \begin{figure}
342     \includegraphics[width=3in]{QuadrupoleFigure}
343     \caption{With the real-space electrostatic methods, the effective
344     quadrupolar tensor, $\mathbf{T}_{\alpha\beta\gamma\delta}(r)$, governing interactions between molecular quadrupoles can be represented by interaction of charges, point dipoles or single point quadrupoles.}
345     \label{fig:quadrupolarFluid}
346     \end{figure}
347     The integral in equation (\ref{gradMaxwell}) can be divided into two parts, $|\textbf{r}-\textbf{r}'|\rightarrow 0 $ and $|\textbf{r}-\textbf{r}'|> 0$. Since the total
348     field gradient due to quadrupolar fluid vanishes at the singularity (see Appendix \ref{singularQuad}), equation (\ref{gradMaxwell}) can be written as,
349     \begin{equation}
350     \partial_\alpha E_\beta(\textbf{r}) = \partial_\alpha {E^o}_\beta(\textbf{r}) +
351     \frac{1}{4\pi \epsilon_o}\int\limits_{|\textbf{r}-\textbf{r}'|> 0 }
352     T_{\alpha\beta\gamma\delta}(|\textbf{r}-\textbf{r}'|)\;{Q}_{\gamma\delta}(\textbf{r}')\;
353     d^3r'.
354     \end{equation}
355     If $\textbf{r} = \textbf{r}'$ is excluded from the integration, the gradient of the electric can be expressed in terms of traceless quadrupole moment as below, \cite{LoganI81}
356     \begin{equation}
357     \partial_\alpha E_\beta(\textbf{r}) = \partial_\alpha {E^o}_\beta(\textbf{r}) + \frac{1}{12\pi \epsilon_o}\int\limits_{|\textbf{r}-\textbf{r}'|> 0 } T_{\alpha\beta\gamma\delta}(|\textbf{r}-\textbf{r}'|)\;{\Theta}_{\gamma\delta}(\textbf{r}')\; d^3r',
358     \end{equation}
359     where $\Theta_{\alpha\beta} = 3Q_{\alpha\beta} - \delta_{\alpha\beta}Tr(Q)$
360     is the traceless quadrupole moment. The total quadrupolar polarization is written as,
361     \begin{equation}
362     {Q}_{\alpha\beta}(\textbf{r}) = \frac{1}{3}\delta_{\alpha\beta}\;Tr({Q})+\epsilon_o {\chi}^*_Q\;\partial_\alpha E_\beta(\textbf{r}),
363     \label{constQaud}
364     \end{equation}
365     In the equation (\ref{constQaud}), $\partial_{\alpha}E_{\beta}$ is Maxwell field gradient and ${\chi}^*_Q$ is the actual quadrupolar susceptibility of the fluid which is different than the proportionality constant $\chi_q $ in the equation (\ref{propConstQuad}). In terms of traceless quadrupole moment, equation (\ref{constQaud}) can be written as,
366     \begin{equation}
367     \frac{1}{3}{\Theta}_{\alpha\beta}(\textbf{r}) = \epsilon_o {\chi}^*_Q \; \partial_\alpha E_\beta (\textbf{r})= \epsilon_o {\chi}^*_Q \left(\partial_\alpha {E^o}_\beta(\textbf{r}) + \frac{1}{12\pi \epsilon_o}\int\limits_{|\textbf{r}-\textbf{r}'|> 0 } T_{\alpha\beta\gamma\delta}(|\textbf{r}-\textbf{r}'|)\;{\Theta}_{\gamma\delta}(\textbf{r}')\; d^3r'\right)
368     \end{equation}
369     For toroidal boundary conditions, both $\partial_\alpha E_\beta$ and
370     $\bar{\Theta}_{\alpha\beta}$ are uniform over the entire space. After
371     performing a Fourier transform (see the Appendix in the Neumann's Paper \cite{Neumann83}) we get,
372     \begin{equation}
373     \frac{1}{3}{{\Theta}}_{\alpha\beta}({\kappa})=
374     \epsilon_o {\chi}^*_Q \;\left[{\partial_\alpha
375     {E^o}_\beta}({\kappa})+ \frac{1}{12\pi
376     \epsilon_o}\;{T}_{\alpha\beta\gamma\delta}({\kappa})\;
377     {{\Theta}}_{\gamma\delta}({\kappa})\right]
378     \end{equation}
379     Since the quadrupolar polarization is in the direction of the applied
380     field, we can write
381     ${{\Theta}}_{\gamma\delta}({\kappa}) =
382     {{\Theta}}_{\alpha\beta}({\kappa})$
383     and
384     ${T}_{\alpha\beta\gamma\delta}({\kappa}) =
385     {T}_{\alpha\beta\alpha\beta}({\kappa})$
386     (note: This is true for orientational polarization but might not be
387     true for polarization due to changes in the atomic distribution. For
388     the orientational case, we only need to work with individual
389     components of the matrix).
390     \begin{equation}
391     \begin{split}
392     \frac{1}{3}{{\Theta}}_{\alpha\beta}({\kappa}) &= \epsilon_o {\chi}^*_Q \left[{\partial_\alpha E^o_\beta}({\kappa})+ \frac{1}{12\pi \epsilon_o}{T}_{\alpha\beta\alpha\beta}({\kappa})\;{{\Theta}}_{\alpha\beta}({\kappa})\right] \\
393     &= \epsilon_o {\chi}^*_Q\;\left(1-\frac{1}{4\pi} {\chi}^*_Q\;
394     {T}_{\alpha\beta\alpha\beta}({\kappa})\right)^{-1}
395     {\partial_\alpha E^o_\beta}({\kappa})
396     \end{split}
397     \label{fourierQuad}
398     \end{equation}
399     Matrix contraction of quadrupole-quadrupole tensor and quadrupole
400     moment gives the identity matrix. Therefore each component can be
401     treated as scalar. If the field gradient is homogeneous over the
402     entire volume, ${\partial_ \alpha E_\beta}({\kappa}) = 0 $ except at
403     $ {\kappa} = 0$, hence it is sufficient to know
404     ${T}_{\alpha\beta\alpha\beta}({\kappa})$ at $ {\kappa} =
405     0$. Therefore equation (\ref{fourierQuad}) can be written as,
406     \begin{equation}
407     \begin{split}
408     \frac{1}{3}{{\Theta}}_{\alpha\beta}({0}) &= \epsilon_o {\chi}^*_Q\; \left(1-\frac{1}{4\pi} {\chi}^*_Q\;{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1} \partial_\alpha E^o_\beta({0})
409     \end{split}
410     \label{fourierQuad2}
411     \end{equation}
412     where $ {T}_{\alpha\beta\alpha\beta}({0})$ can be evaluated as,
413     \begin{equation}
414     {T}_{\alpha\beta\alpha\beta}({0}) = \int {T}_{\alpha\beta\alpha\beta}\;(\textbf{r})d^3r
415     \label{realTensorQaud}
416     \end{equation}
417    
418     In terms of traced quadrupole moment equation (\ref{fourierQuad2}) can be written as,
419     \begin{equation}
420     {{Q}}_{\alpha\beta} = \frac{1}{3}\delta_{\alpha\beta}\;Tr({Q}) + \epsilon_o\; {\chi}^*_Q\left(1-\frac{1}{4\pi} {\chi}^*_Q\;{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1}\; \partial_\alpha E^o_\beta
421     \label{tracedConstQuad}
422     \end{equation}
423     Comparing (\ref{tracedConstQuad}) and (\ref{flucQuad}) we get,
424     \begin{equation}
425     \begin{split}
426     &\frac{\braket{{Q^2}} - \braket{Q}^2}{15 \epsilon_o Vk_BT}\; =\; {\chi}^*_Q\;\left(1-\frac{1}{4\pi} {\chi}^*_Q\;{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1}, \\
427     &{\chi}^*_Q \;=\; \frac{\braket{{Q^2}} - \braket{Q}^2}{15 \epsilon_o Vk_BT}\left(1 + \frac{1}{4\pi} \frac{\braket{{Q^2}} - \braket{Q}^2}{15 \epsilon_o Vk_BT}\;{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1}
428     \end{split}
429     \end{equation}
430     Finally the quadrupolar susceptibility cab be written in terms of quadrupolar correction factor ($A_{quad}$) as below,
431     \begin{equation}
432     {\chi}^*_Q \;=\; \frac{\braket{{Q^2}} - \braket{Q}^2}{15 \epsilon_o Vk_BT}\left(1 + \frac{\braket{{Q^2}} - \braket{Q}^2}{15 \epsilon_o Vk_BT}\; A_{quad}\right)^{-1}
433     \end{equation}
434     where $A_{quad} = \frac{1}{4\pi}\int {T}_{\alpha\beta\alpha\beta}\;(\textbf{r})d^3r $ has dimension of the $length^{-2}$ is different for different cutoff methods which is listed in Table \ref{tab:B}. The dielectric constant associated with the quadrupolar susceptibility is defined as,\cite{Ernst92}
435    
436     \begin{equation}
437     \epsilon = 1 + \chi^*_Q\; G,
438     \end{equation}
439     where $G = \frac{\displaystyle\int_V |\partial_\alpha E^o_\beta|^2 d^3r}{\displaystyle\int_V{|E^o|}^2 d^3r}$ is a geometrical factor depends on the nature of the external field perturbation. This is true when the quadrupolar fluid is homogeneous over the sample. Since quadrupolar molecule couple with the gradient of the field, the distribution of the quadrupoles is inhomogeneous for varying field gradient. Hence the distribution function should also be taken into account to calculate actual geometrical factor in the presence of non-uniform gradient field. Therefore,
440     \begin{equation}
441     G = \frac{\displaystyle\int_V\; g(r, \theta, \phi)\; |\partial_\alpha E^o_\beta|^2 d^3r}{\displaystyle\int_V{|E^o|}^2 d^3r}
442     \end{equation}
443     where $g(r,\theta, \phi)$ is a distribution function of the quadrupoles in with respect to origin at the center of line joining two probe charges.
444     \begin{table}
445     \caption{Expressions for the quadrupolar correction factor ($A$) for the real-space electrostatic methods in terms of the damping parameter
446     ($\alpha$) and the cutoff radius ($r_c$). The dimension of the correction factor is $ length^{-2}$ in case of quadrupolar fluid.}
447     \label{tab:B}
448     \centering
449     \resizebox{\columnwidth}{!}{%
450    
451     \begin{tabular}{l|c|c|c|c|}
452    
453     Method & $A_\mathrm{charges}$ & $A_\mathrm{dipoles}$ &$A_\mathrm{quadrupoles}$ \\\hline
454     Spherical Cutoff (SC) & $ -\frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $ -\frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $ -\frac{8 {\alpha}^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ \\
455     Shifted Potental (SP) & $ -\frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $- \frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$& $ -\frac{16 \alpha^7 {r_c}^5}{9\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ \\
456     Gradient-shifted (GSF) & $- \frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & 0 & $-\frac{4{\alpha}^7{r_c}^5 }{9\sqrt{\pi}}e^{-\alpha^2 r_c^2}(-1+2\alpha ^2 r_c^2)$\\
457     Taylor-shifted (TSF) & $ -\frac{8 \alpha^5 {r_c}^3}{3\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $\frac{4\;\mathrm{erfc(\alpha r_c)}}{{r_c}^2} + \frac{8 \alpha}{3\sqrt{\pi}r_c}e^{-\alpha^2 {r_c}^2}\left(3+ 2 \alpha^2 {r_c}^2 + \alpha^4 {r_c}^4\right) $ & $\frac{10\;\mathrm{erfc}(\alpha r_c )}{{r_c}^2} + \frac{4{\alpha}}{9\sqrt{\pi}{r_c}}e^{-\alpha^2 r_c^2}\left(45 + 30\alpha ^2 {r_c}^2 + 12\alpha^4 {r_c}^4 + 3\alpha^6 {r_c}^6 + 2 \alpha^8 {r_c}^8\right)$ \\\hline
458     \end{tabular}%
459     }
460     \end{table}
461     \section{Methodology}
462    
463     \section{Result}
464    
465     \section{Conclusion}
466    
467     \newpage
468    
469     \appendix
470     \section{Point-multipolar interactions with a spatially-varying electric field}
471    
472     We can treat objects $a$, $b$, and $c$ containing embedded collections
473     of charges. When we define the primitive moments, we sum over that
474     collections of charges using a local coordinate system within each
475     object. The point charge, dipole, and quadrupole for object $a$ are
476     given by $C_a$, $\mathbf{D}_a$, and $\mathsf{Q}_a$, respectively.
477     These are the primitive multipoles which can be expressed as a
478     distribution of charges,
479     \begin{align}
480     C_a =&\sum_{k \, \text{in }a} q_k , \label{eq:charge} \\
481     D_{a\alpha} =&\sum_{k \, \text{in }a} q_k r_{k\alpha}, \label{eq:dipole}\\
482     Q_{a\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in } a} q_k
483     r_{k\alpha} r_{k\beta} . \label{eq:quadrupole}
484     \end{align}
485     Note that the definition of the primitive quadrupole here differs from
486     the standard traceless form, and contains an additional Taylor-series
487     based factor of $1/2$. In Paper 1, we derived the forces and torques
488     each object exerts on the others.
489    
490     Here we must also consider an external electric field that varies in
491     space: $\mathbf E(\mathbf r)$. Each of the local charges $q_k$ in
492     object $a$ will then experience a slightly different field. This
493     electric field can be expanded in a Taylor series around the local
494     origin of each object. A different Taylor series expansion is carried
495     out for each object.
496    
497     For a particular charge $q_k$, the electric field at that site's
498     position is given by:
499     \begin{equation}
500     E_\gamma + \nabla_\delta E_\gamma r_{k \delta}
501     + \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta}
502     r_{k \varepsilon} + ...
503     \end{equation}
504     Note that the electric field is always evaluated at the origin of the
505     objects, and treating each object using point multipoles simplifies
506     this greatly.
507    
508     To find the force exerted on object $a$ by the electric field, one
509     takes the electric field expression, and multiplies it by $q_k$, and
510     then sum over all charges in $a$:
511    
512     \begin{align}
513     F_\gamma &= \sum_{k \textrm{~in~} a} q_k \lbrace E_\gamma + \nabla_\delta E_\gamma r_{k \delta}
514     + \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta}
515     r_{k \varepsilon} + ... \rbrace \\
516     &= C_a E_\gamma + D_{a \delta} \nabla_\delta E_\gamma
517     + Q_{a \delta \varepsilon} \nabla_\delta \nabla_\varepsilon E_\gamma +
518     ...
519     \end{align}
520    
521     Similarly, the torque exerted by the field on $a$ can be expressed as
522     \begin{align}
523     \tau_\alpha &= \sum_{k \textrm{~in~} a} (\mathbf r_k \times q_k \mathbf E)_\alpha \\
524     & = \sum_{k \textrm{~in~} a} \epsilon_{\alpha \beta \gamma} q_k
525     r_{k\beta} E_\gamma(\mathbf r_k) \\
526     & = \epsilon_{\alpha \beta \gamma} D_\beta E_\gamma
527     + 2 \epsilon_{\alpha \beta \gamma} Q_{\beta \delta} \nabla_\delta
528     E_\gamma + ...
529     \end{align}
530    
531     The last term is essentially identical with form derived by Torres del
532     Castillo and M\'{e}ndez Garrido,\cite{Torres-del-Castillo:2006uo} although their derivation
533     utilized a traceless form of the quadrupole that is different than the
534     primitive definition in use here. We note that the Levi-Civita symbol
535     can be eliminated by utilizing the matrix cross product in an
536     identical form as in Ref. \onlinecite{Smith98}:
537     \begin{equation}
538     \left[\mathsf{A} \times \mathsf{B}\right]_\alpha = \sum_\beta
539     \left[\mathsf{A}_{\alpha+1,\beta} \mathsf{B}_{\alpha+2,\beta}
540     -\mathsf{A}_{\alpha+2,\beta} \mathsf{B}_{\alpha+1,\beta}
541     \right]
542     \label{eq:matrixCross}
543     \end{equation}
544     where $\alpha+1$ and $\alpha+2$ are regarded as cyclic permuations of
545     the matrix indices. In table \ref{tab:UFT} we give compact
546     expressions for how the multipole sites interact with an external
547     field that has exhibits spatial variations.
548    
549     \begin{table}
550     \caption{Potential energy $(U)$, force $(\mathbf{F})$, and torque
551     $(\mathbf{\tau})$ expressions for a multipolar site embedded in an
552     electric field with spatial variations, $\mathbf{E}(\mathbf{r})$.
553     \label{tab:UFT}}
554     \begin{tabular}{r|ccc}
555     & Charge & Dipole & Quadrupole \\ \hline
556     $U$ & $C \phi(\mathbf{r})$ & $-\mathbf{D} \cdot \mathbf{E}(\mathbf{r})$ & $- \mathsf{Q}:\nabla \mathbf{E}(\mathbf{r})$ \\
557     $\mathbf{F}$ & $C \mathbf{E}(\mathbf{r})$ & $+\mathbf{D} \cdot \nabla \mathbf{E}(\mathbf{r})$ & $+\mathsf{Q} : \nabla\nabla\mathbf{E}(\mathbf{r})$ \\
558     $\mathbf{\tau}$ & & $\mathbf{D} \times \mathbf{E}(\mathbf{r})$ & $+2 \mathsf{Q} \times \nabla \mathbf{E}(\mathbf{r})$
559     \end{tabular}
560     \end{table}
561     \section{Gradient of the field due to quadrupolar polarization}
562     \label{singularQuad}
563     In this section, we will discuss the gradient of the field produced by
564     quadrupolar polarization. For this purpose, we consider a distribution
565     of charge ${\rho}(r)$ which gives rise to an electric field
566     $\vec{E}(r)$ and gradient of the field $\vec{\nabla} \vec{E}(r)$
567     throughout space. The total gradient of the electric field over volume
568     due to the all charges within the sphere of radius $R$ is given by
569     (cf. Jackson equation 4.14):
570     \begin{equation}
571     \int_{r<R} \vec{\nabla}\vec{E}\;d^3r = -\int_{r=R} R^2 \vec{E}\;\hat{n}\; d\Omega
572     \label{eq:8}
573     \end{equation}
574     where $d\Omega$ is the solid angle and $\hat{n}$ is the normal vector
575     of the surface of the sphere which is equal to
576     $sin[\theta]cos[\phi]\hat{x} + sin[\theta]sin[\phi]\hat{y} +
577     cos[\theta]\hat{z}$
578     in spherical coordinates. For the charge density ${\rho}(r')$, the
579     total gradient of the electric field can be written as (cf. Jackson
580     equation 4.16),
581     \begin{equation}
582     \int_{r<R} \vec{\nabla}\vec{E}\; d^3r=-\int_{r=R} R^2\; \vec{\nabla}\Phi\; \hat{n}\; d\Omega =-\frac{1}{4\pi\;\epsilon_o}\int_{r=R} R^2\; \vec{\nabla}\;\left(\int \frac{\rho(r')}{|\vec{r}-\vec{r'}|}\;d^3r'\right) \hat{n}\; d\Omega
583     \label{eq:9}
584     \end{equation}
585     The radial function in the equation (\ref{eq:9}) can be expressed in
586     terms of spherical harmonics as (cf. Jackson equation 3.70),
587     \begin{equation}
588     \frac{1}{|\vec{r} - \vec{r'}|} = 4\pi \sum_{l=0}^{\infty}\sum_{m=-l}^{m=l}\frac{1}{2l+1}\;\frac{{r^l_<}}{{r^{l+1}_>}}\;{Y^*}_{lm}(\theta', \phi')\;Y_{lm}(\theta, \phi)
589     \label{eq:10}
590     \end{equation}
591     If the sphere completely encloses the charge density then $ r_< = r'$ and $r_> = R$. Substituting equation (\ref{eq:10}) into (\ref{eq:9}) we get,
592     \begin{equation}
593     \begin{split}
594     \int_{r<R} \vec{\nabla}\vec{E}\;d^3r &=-\frac{R^2}{\epsilon_o}\int_{r=R} \; \vec{\nabla}\;\left(\int \rho(r')\sum_{l=0}^{\infty}\sum_{m=-l}^{m=l}\frac{1}{2l+1}\;\frac{{r'^l}}{{R^{l+1}}}\;{Y^*}_{lm}(\theta', \phi')\;Y_{lm}(\theta, \phi)\;d^3r'\right) \hat{n}\; d\Omega \\
595     &= -\frac{R^2}{\epsilon_o}\sum_{l=0}^{\infty}\sum_{m=-l}^{m=l}\frac{1}{2l+1}\;\int \rho(r')\;{r'^l}\;{Y^*}_{lm}(\theta', \phi')\left(\int_{r=R}\vec{\nabla}\left({R^{-(l+1)}}\;Y_{lm}(\theta, \phi)\right)\hat{n}\; d\Omega \right)d^3r
596     '
597     \end{split}
598     \label{eq:11}
599     \end{equation}
600     The gradient of the product of radial function and spherical harmonics
601     is given by (cf. Arfken, p.811 eq. 16.94):
602     \begin{equation}
603     \begin{split}
604     \vec{\nabla}\left[ f(r)\;Y_{lm}(\theta, \phi)\right] = &-\left(\frac{l+1}{2l+1}\right)^{1/2}\; \left[\frac{\partial}{\partial r}-\frac{l}{r} \right]f(r)\; Y_{l, l+1, m}(\theta, \phi)\\ &+ \left(\frac{l}{2l+1}\right)^{1/2}\left[\frac
605     {\partial}{\partial r}+\frac{l}{r} \right]f(r)\; Y_{l, l-1, m}(\theta, \phi).
606     \end{split}
607     \label{eq:12}
608     \end{equation}
609     Using equation (\ref{eq:12}) we get,
610     \begin{equation}
611     \vec{\nabla}\left({R^{-(l+1)}}\;Y_{lm}(\theta, \phi)\right) = [(l+1)(2l+1)]^{1/2}\; Y_{l,l+1,m}(\theta, \phi) \; \frac{1}{R^{l+2}},
612     \label{eq:13}
613     \end{equation}
614     where $ Y_{l,l+1,m}(\theta, \phi)$ is the vector spherical harmonics
615     which can be expressed in terms of spherical harmonics as shown in
616     below (cf. Arfkan p.811),
617     \begin{equation}
618     Y_{l,l+1,m}(\theta, \phi) = \sum_{m_1, m_2} C(l+1,1,l|m_1,m_2,m)\; {Y_{l+1}}^{m_1}(\theta,\phi)\; \hat{e}_{m_2},
619     \label{eq:14}
620     \end{equation}
621     where $C(l+1,1,l|m_1,m_2,m)$ is a Clebsch-Gordan coefficient and
622     $\hat{e}_{m_2}$ is a spherical tensor of rank 1 which can be expressed
623     in terms of Cartesian coordinates,
624     \begin{equation}
625     {\hat{e}}_{+1} = - \frac{\hat{x}+i\hat{y}}{\sqrt{2}},\quad {\hat{e}}_{0} = \hat{z},\quad and \quad {\hat{e}}_{-1} = \frac{\hat{x}-i\hat{y}}{\sqrt{2}}
626     \label{eq:15}
627     \end{equation}
628     The normal vector $\hat{n} $ can be expressed in terms of spherical tensor of rank 1 as shown in below,
629     \begin{equation}
630     \hat{n} = \sqrt{\frac{4\pi}{3}}\left(-{Y_1}^{-1}{\hat{e}}_1 -{Y_1}^{1}{\hat{e}}_{-1} + {Y_1}^{0}{\hat{e}}_0 \right)
631     \label{eq:16}
632     \end{equation}
633     The surface integral of the product of $\hat{n}$ and
634     ${Y_{l+1}}^{m_1}(\theta, \phi)$ gives,
635     \begin{equation}
636     \begin{split}
637     \int \hat{n}\;{Y_{l+1}}^{m_1}\;d\Omega &= \int \sqrt{\frac{4\pi}{3}}\left(-{Y_1}^{-1}{\hat{e}}_1 -{Y_1}^{1}{\hat{e}}_{-1} + {Y_1}^{0}{\hat{e}}_0 \right)\;{Y_{l+1}}^{m_1}\; d\Omega \\
638     &= \int \sqrt{\frac{4\pi}{3}}\left({{Y_1}^{1}}^* {\hat{e}}_1 +{{Y_1}^{-1}}^* {\hat{e}}_{-1} + {{Y_1}^{0}}^* {\hat{e}}_0 \right)\;{Y_{l+1}}^{m_1}\; d\Omega \\
639     &= \sqrt{\frac{4\pi}{3}}\left({\delta}_{l+1, 1}\;{\delta}_{1, m_1}\;{\hat{e}}_1 + {\delta}_{l+1, 1}\;{\delta}_{-1, m_1}\;{\hat{e}}_{-1}+ {\delta}_{l+1, 1}\;{\delta}_{0, m_1} \;{\hat{e}}_0\right),
640     \end{split}
641     \label{eq:17}
642     \end{equation}
643     where ${Y_{l}}^{-m} = (-1)^m\;{{Y_{l}}^{m}}^* $ and
644     $ \int {{Y_{l}}^{m}}^*\;{Y_{l'}}^{m'}\;d\Omega =
645     \delta_{ll'}\delta_{mm'} $.
646     Non-vanishing values of equation \ref{eq:17} require $l = 0$,
647     therefore the value of $ m = 0 $. Since the values of $ m_1$ are -1,
648     1, and 0 then $m_2$ takes the values 1, -1, and 0, respectively
649     provided that $m = m_1 + m_2$. Equation \ref{eq:11} can therefore be
650     modified,
651     \begin{equation}
652     \begin{split}
653     \int_{r<R} \vec{\nabla}\vec{E}\;d^3r = &- \sqrt{\frac{4\pi}{{3}}}\;\frac{1}{\epsilon_o}\int \rho(r')\;{Y^*}_{00}(\theta', \phi')[ C(1, 1, 0|-1,1,0)\;{\hat{e}_{-1}}{\hat{e}_{1}}\\ &+ C(1, 1, 0|-1,1,0)\;{\hat{e}_{1}}{\hat{e}_{-1}}+C(
654     1, 1, 0|0,0,0)\;{\hat{e}_{0}}{\hat{e}_{0}} ]\; d^3r'.
655     \end{split}
656     \label{eq:18}
657     \end{equation}
658     After substituting ${Y^*}_{00} = \frac{1}{\sqrt{4\pi}} $ and using the
659     values of the Clebsch-Gorden coefficients: $ C(1, 1, 0|-1,1,0) =
660     \frac{1}{\sqrt{3}}, \; C(1, 1, 0|-1,1,0)= \frac{1}{\sqrt{3}}$ and $
661     C(1, 1, 0|0,0,0) = -\frac{1}{\sqrt{3}}$ in equation \ref{eq:18} we
662     obtain,
663     \begin{equation}
664     \begin{split}
665     \int_{r<R} \vec{\nabla}\vec{E}\;d^3r &= -\sqrt{\frac{4\pi}{{3}}}\;\frac{1}{\epsilon_o}\int \rho(r')\;d^3r'\left({\hat{e}_{-1}}{\hat{e}_{1}}+{\hat{e}_{1}}{\hat{e}_{-1}}-{\hat{e}_{0}}{\hat{e}_{0}}\right)\\
666     &= - \sqrt{\frac{4\pi}{{3}}}\;\frac{1}{\epsilon_o}\;C_{total}\;\left({\hat{e}_{-1}}{\hat{e}_{1}}+{\hat{e}_{1}}{\hat{e}_{-1}}-{\hat{e}_{0}}{\hat{e}_{0}}\right).
667     \end{split}
668     \label{eq:19}
669     \end{equation}
670     Equation (\ref{eq:19}) gives the total gradient of the field over a
671     sphere due to the distribution of the charges. For quadrupolar fluids
672     the total charge within a sphere is zero, therefore
673     $ \int_{r<R} \vec{\nabla}\vec{E}\;d^3r = 0 $. Hence the quadrupolar
674     polarization produces zero net gradient of the field inside the
675     sphere.
676    
677    
678    
679     \newpage
680    
681     \bibliography{multipole}
682    
683     \end{document}