trunk/multipole/dielectric.tex

\documentclass[%
 aip,
 jcp,
 amsmath,amssymb,
preprint,%
% reprint,%
%author-year,%
%author-numerical,%
]{revtex4-1}

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{multirow}
\usepackage{bm}% bold math
\usepackage{natbib}
\usepackage{times}
\usepackage{mathptmx}
\usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions
\usepackage{url}
\usepackage{braket}

\begin{document}

\title{Real space electrostatics for multipoles. III. Dielectric Properties}

\author{Madan Lamichhane}
\affiliation{Department of Physics, University
of Notre Dame, Notre Dame, IN 46556}
\author{Thomas Parsons}
\affiliation{Department of Chemistry and Biochemistry, University
of Notre Dame, Notre Dame, IN 46556}
\author{Kathie E. Newman}
\affiliation{Department of Physics, University
of Notre Dame, Notre Dame, IN 46556}
\author{J. Daniel Gezelter}
\email{gezelter@nd.edu.}
\affiliation{Department of Chemistry and Biochemistry, University
of Notre Dame, Notre Dame, IN 46556}

\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\begin{abstract}
  Note: This manuscript is a work in progress.   

  We report on the dielectric properties of the shifted potential
  (SP), gradient shifted force (GSF), and Taylor shifted force (TSF)
  real-space methods for multipole interactions that were developed in
  the first two papers in this series.  We find that some subtlety is
  required for computing dielectric properties with the real-space
  methods, particularly when using the common fluctuation formulae.
  Three distinct methods for computing the dielectric constant are
  investigated, including the standard fluctuation formulae,
  potentials of mean force between solvated ions, and direct
  measurement of linear solvent polarization in response to applied
  fields and field gradients.
\end{abstract}

\maketitle

\section{Introduction}

Over the past several years, there has been increasing interest in
pairwise methods for correcting electrostatic interactions in computer
simulations of condensed molecular
systems.\cite{Wolf99,Zahn02,Kast03,Beckd.A.C._Bi0486381,Ma05,Fennell06}
These techniques were developed from the observations and efforts of
Wolf {\it et al.}  and their work towards an $\mathcal{O}(N)$
Coulombic sum.\cite{Wolf99} Wolf's method of cutoff neutralization is
able to obtain excellent agreement with Madelung energies in ionic
crystals.\cite{Wolf99} 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.  Of particular interest is the
static dielectric constant, $\epsilon$.  Using the Ewald sum under
tin-foil boundary conditions, $\epsilon$ can be calculated for systems
of non-polarizable substances via
\begin{equation}
\epsilon = 1 + \frac{\langle M^2\rangle}{3\epsilon_0k_\textrm{B}TV},
\label{eq:staticDielectric}
\end{equation}
where $\epsilon_0$ is the permittivity of free space and $\langle
M^2\rangle$ is the fluctuation of the system dipole
moment.\cite{Allen:1989fp} The numerator in the fractional term in equation
(\ref{eq:staticDielectric}) is the fluctuation of the simulation-box
dipole moment, identical to the quantity calculated in the
finite-system Kirkwood $g$ factor ($G_k$):
\begin{equation}
G_k = \frac{\langle M^2\rangle}{N\mu^2},
\label{eq:KirkwoodFactor}
\end{equation}
where $\mu$ is the dipole moment of a single molecule of the
homogeneous system.\cite{Neumann:1983mz,Neumann:1983yq,Neumann84,Neumann85} The
fluctuation term in both equation (\ref{eq:staticDielectric}) and
(\ref{eq:KirkwoodFactor}) is calculated as follows,
\begin{equation}
\begin{split}
\langle M^2\rangle &= \langle\bm{M}\cdot\bm{M}\rangle 
                        - \langle\bm{M}\rangle\cdot\langle\bm{M}\rangle \\
                   &= \langle M_x^2+M_y^2+M_z^2\rangle
                        - (\langle M_x\rangle^2 + \langle M_y\rangle^2 
                                + \langle M_z\rangle^2).        
\end{split}
\label{eq:fluctBoxDipole}
\end{equation}
This fluctuation term can be accumulated during a simulation; however,
it converges rather slowly, thus requiring multi-nanosecond simulation
times.\cite{Horn04} In the case of tin-foil boundary conditions, the
dielectric/surface term of the Ewald summation is equal to zero.

\section{Dielectric constants for real-space electrostatics}

The new real-space methods require some care when computing dielectric
properties because the interaction between molecular dipoles depends
on the level of approximation being utilized. 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}).
\begin{figure}
\includegraphics[width=3in]{DielectricFigure}
\caption{With the real-space electrostatic methods, the effective
  dipole tensor, $\mathbf{T}$, governing interactions between
  molecular dipoles is not the same for charge-charge interactions as
  for point dipoles.}
\label{fig:stockmayer}
\end{figure}
In the case where point charges are interacting via an electrostatic
kernel, $v(r)$, the effective {\it molecular} dipole tensor,
$\mathbf{T}$ is obtained from two successive applications of the
gradient operator to the electrostatic kernel,
\begin{equation}
\mathbf{T}_{\alpha \beta}(r) =  \nabla_\alpha \nabla_\beta \left(v(r)\right) = \delta_{\alpha \beta}
\left(\frac{1}{r} v^\prime(r) \right) + \frac{r_{\alpha}
  r_{\beta}}{r^2} \left( v^{\prime \prime}(r) - \frac{1}{r}
  v^{\prime}(r) \right)
\end{equation}
where $v(r)$ may be either the bare kernel ($1/r$) or one of the
modified (Wolf or DSF) kernels.  This tensor describes the effective
interaction between molecular dipoles ($\mathbf{D}$) in Gaussian
units as $-\mathbf{D} \cdot \mathbf{T} \cdot \mathbf{D}$.

When utilizing the new real-space methods for point dipoles, the
tensor is explicitly constructed,
\begin{equation}
\mathbf{T}_{\alpha \beta}(r)  =  \delta_{\alpha \beta} v_{21}(r) +
\frac{r_{\alpha} r_{\beta}}{r^2} v_{22}(r) 
\end{equation}
where the functions $v_{21}(r)$ and $v_{22}(r)$ depend on the level of
the approximation. Although the Taylor-shifted (TSF) and
gradient-shifted (GSF) models produce to the same $v(r)$ function for
point charges, they have distinct forms for the dipole-dipole
interactions.

Neumann and Steinhauser showed~\cite{Neumann:1983mz,Neumann:1983yq}
that the relative dielectric permittivity $\epsilon$ is given by the
general fluctuation formula,
\begin{equation}
\frac{\epsilon - 1}{\epsilon +2} = \frac{4 \pi}{3} \frac{\left<M^2\right>}{3 V k_B
  T} \left( 1 - \frac{3}{4 \pi} \frac{\epsilon -1}{\epsilon + 2}
  \tilde{T}(0) \right)
\end{equation}
where $\left<M^2\right>$ is the mean fluctuation of the square of the
net dipole moment of the simulation cell, $V$ is the volume of the
cell, and $\tilde{T}(0) = \tilde{T̃}_{xx}(0) = \tilde{T̃}_{yy}(0) =
\tilde{T̃}_{zz}(0)$ is the $\mathbf{k} = 0$ limit of Fourier transform
of the diagonal term of the effective molecular dipolar tensor,
\begin{equation}
\tilde{T}(0) = \int_V \mathbf{T}(\mathbf{r})  d\mathbf{r}
\end{equation}
where the integral is carried out over the relevant geometry for the
interaction.  For the real-space methods, the integration can be done
easily up to the imposed cutoff radius, while for the Ewald sum, the
results also depend on details of the $k$-space
calculation.\cite{Neumann:1983yq}

For molecules composed of point charges interacting via the bare
$(1/r)$ kernel, the simple cutoff (SC) and minimum image (MI)
approaches give $\tilde{T}(0) = 0$ which means that the
Clausius-Mosotti equation,
\begin{equation}
\frac{4\pi}{3}\frac{\left< M^2 \right>}{3 V k_B T} = \frac{\epsilon
  -1}{\epsilon+2} 
\end{equation}
is the relevant fluctuation formula for the relative permittivity.

Zahn {\it et al.}\cite{Zahn02} showed that for the damped shifted
force charge-charge kernel, $\tilde{T}(0) = 4 \pi / 3$. This was later
generalized by Izvekov {\it et al.} for all point-charge kernels which
have forces that go to zero at a cutoff radius and which maintain a
pole of first order at $r=0$.\cite{Izvekov:2008wo} When $\tilde{T}(0)
= 4 \pi/ 3$, the expression for the dielectric constant reduces to the
widely-used {\it conducting boundary} formula,
\begin{equation}
\frac{4\pi}{3}\frac{\left< M^2 \right>}{3 V k_B T} = \frac{\epsilon -
  1}{3}
\end{equation}
It is convenient to define a quantity $A = \frac{3}{4 \pi}
\tilde{T}(0)$, and to combine all of the fluctuation formulae
together:
\begin{equation}
\frac{4 \pi}{3} \frac{\left< M^2 \right>}{3 V k_B T} =
\left\{\begin{array}{ll}
\frac{\epsilon-1}{\epsilon+2} \left[1- \frac{\epsilon-1}{\epsilon+2}
  A\right]^{-1} & \mathrm{General~Case} \\
\frac{\epsilon-1}{\epsilon+2} & A \rightarrow 0 \mathrm{~limit} \\
\frac{\epsilon-1}{3} & A \rightarrow 1 \mathrm{~limit}
\end{array}\right.
\end{equation}

The Clausius-Mossotti $(A\rightarrow 0)$ approach is subject to noise
and error magnification,\cite{Allen:1989fp} and the conducting
boundary approach $(A \rightarrow 1)$ has become widely used as a
result.  If the electrostatic method being used has $A < 1$, the
relative dielectric permittivity $\epsilon$ can be estimated once the
conducting boundary fluctuation result has been found,
\begin{equation}
\epsilon = \frac{(A+2) (\epsilon_\text{CB}-1)+3} {(A-1) (\epsilon_\text{CB}-1)+3}
\end{equation}
Note that this expression becomes quite numerically sensitive when the
value of $A$ deviates significantly from 1.

We have derived a set of expressions for the value of $A$ for the new
real space methods that are shown in Table \ref{tab:A}.

\begin{table}
  \caption{Expressions for the dielectric correction factor ($A$) for the
    real-space electrostatic methods in terms of the damping parameter
    ($\alpha$) and the cutoff radius ($r_c$).  The Ewald-Kornfeld result 
    derived in Refs. \onlinecite{Adams:1980rt,Adams:1981fr,Neumann:1983yq} is shown for comparison using the Ewald 
    convergence parameter ($\kappa$) and the real-space cutoff value ($r_c$). }
\label{tab:A}
\resizebox{\columnwidth}{!}{%
\begin{tabular}{l|c|c|c|c|c|}
       & \multicolumn{3}{c|}{damped} & \multicolumn{2}{c|}{undamped} \\ \cline{2-4} \cline{5-6}
Method & $A_\mathrm{charges}$ & $A_\mathrm{dipoles}$ & $A_\mathrm{quad} (length^{-2})$ & $A_\mathrm{charges}$ & $A_\mathrm{dipoles}$ \\
\hline
Spherical Cutoff (SC) & $\mathrm{erf}(r_c \alpha) - \frac{2 \alpha r_c}{\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & 0 &$-\frac{8{\alpha}^5{r_c}^3 }{5\sqrt{\pi}}e^{-{\alpha}^2{r_c}^2 }$  &0 &0\\
Shifted Potental (SP) & $ \mathrm{erf}(r_c \alpha) - \frac{2 \alpha r_c}{\sqrt{\pi}} e^{-\alpha^2 r_c^2}$ & $\mathrm{erf}(\alpha r_c)-\frac{2 \alpha  r_c }{\sqrt{\pi }} \left(1 + \frac{2 \alpha^2 r_c^2}{3} \right) e^{-\alpha ^2 r_c^2} $&$-\frac{16{\alpha}^7{r_c}^5 }{15\sqrt{\pi}}e^{-{\alpha}^2{r_c}^2 }$  &0 &0\\
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} $ &$-\frac{4{\alpha}^7{r_c}^5 }{15\sqrt{\pi}}(-1+2\alpha ^2 r_c^2)e^{-\alpha^2 r_c^2}$ &1 &0\\
Taylor-shifted  (TSF) & 1 & 1 &$\frac{6\;\mathrm{erf}(r_c \alpha)}{{r_c}^2} + \frac{4{\alpha}{r_c}}{15\sqrt{\pi}{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)e^{-\alpha^2 r_c^2}$ &1 & 1\\ 
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
\end{tabular}%
}
\end{table}

\section{Boltzman average of the dipole polarization density}
Consider a dipolar fluid with polarization density $\mathbf{P}(r)$ at
equlibrium at temperature $T$. The average of polarization density is
given by,
\begin{equation}
\left<\mathbf{P}(r)\right> = \frac{\int\mathbf{P}(r) e^{-H/k_B T} d
  \Gamma^N}{\int e^{-H/k_B T} d \Gamma^N}
\end{equation}
where $d \Gamma^N = dr_1...dr_N dp_1...dp_N$ is the volume element in
phase space, and $H$ is the Hamiltonian of the $N$-particle system in
the absence of applied external field.  The Hamiltonian of the system
is given by,
\begin{equation}
H = \frac{1}{2}\sum_{i=1}^N m_i {\dot{\mathbf{r}}_i}^2 + \sum_{i<j} U(\mathbf{r}_i - \mathbf{r}_j)
\end{equation}
where $m_i$ is mass of site $i$, and $U(\mathbf{r}_i - \mathbf{r}_j)$
is the interaction potential between sites $i$ and $j$.

Let $E'$ be the applied external electric field.  The Hamiltonian for
the polarized system then becomes
\begin{equation}
H' = H - {M}_\alpha {E'}_\alpha
\end{equation}
where $\mathbf{M}$ is the total dipole moment in zero field. The
polarization density changes in the presence of the external field
$\mathbf{E'}$. The modified polarization is given by
\begin{equation}
\left< P'_\alpha(r) \right>_{E'} = \frac{\int P_\alpha(r) e^{-H'/k_B
    T} d \Gamma^N}{\int e^{-H'/k_B T} d \Gamma^N} \equiv \frac{F}{G}
\end{equation}
For small values of the applied electric field we can expand the
polarization density as,
\begin{equation}
\left< P'_\alpha(r) \right>_{E'} = \left< P'_\alpha(r) \right> + \frac{\partial
  \left< P'_\alpha(r) \right>}{\partial E'_\beta }\bigg|_{E'=0} E'_\beta
\end{equation}
and 
\begin{equation}
\frac{\partial{\left<P'_\alpha(r)\right>}}{\partial{E'_\beta}}\bigg|_{E'=0} = \frac{F'G-FG'}{G^2}
\end{equation}
where
\begin{equation}
\begin{split}
F' & = \frac{1}{k_BT} \int {P}_\alpha M_\beta\; e^{-H/k_B T} d\Gamma^N\\
\textit{and}  \\
G' & = \frac{1}{k_BT} \int M_\beta\; e^{-H/k_B T} d\Gamma^N
\end{split}
\end{equation}
Therefore 
\begin{equation}
\braket{{P'_\alpha(r)}}_{E'}=\braket{{P'_\alpha(r)}}+\frac{1}{k_B T}(\braket{P_{\alpha}M_{\beta}}-\braket{P_{\alpha}}\braket{M_{\beta}}){E'}_\beta,
\end{equation}
where $\braket{A} = \frac{\int A\; e^{-H/k_B T}d\Gamma^N}{\int e^{-H/k_B T}d\Gamma^N}$ is the ensemble average of a physical quantity $A$ at temperature $T$. Now the polarization density due to applied external field is given by,  
\begin{equation}
 {\Delta P}_\alpha =\frac{1}{k_B T}\left(\braket{P_{\alpha}M_{\beta}}-\braket{P_{\alpha}}\braket{M_{\beta}}\right){E'}_\beta.
\end{equation}
Therefore the total polarization for the system can be obtained by integrating polarization density over the volume,
\begin{equation}
 {\Delta M}_\alpha =\frac{1}{k_B T} \left(\braket{M_{\alpha}M_{\beta}}-\braket{M_{\alpha}}\braket{M_{\beta}}\right){E'}_\beta.
\end{equation}
If we consider both applied field and polarization in the z-direction,
\begin{equation}
 {\Delta M}_z =\frac{1}{k_B T}\left(\braket{{M_{z}}^2}-\braket{M_{z}}^2\right){E'}_z.
\end{equation}
In general applied field $ \braket{{M_{z}}^2} = \frac{1}{3}\braket{{\mathbf{M}}^2}$ and $ \braket{{M_{z}}} = \frac{1}{3}\braket{{\mathbf{M}}}$ therefore,
\begin{equation}
{\Delta P}_z = \frac{{\Delta M}_z}{V} = \epsilon_o \frac{\braket{{\mathbf{M}}^2}-\braket{{\mathbf{M}}}^2}{3 \epsilon_o \; k_B T V}{E'}_z.
\label{eq:fluc}
\end{equation}
Comparing equation (\ref{eq:fluc}) with $ \Delta P_z = \epsilon_o \chi \ {E'}_z$ the fluctuation formula for susceptibility ($\chi$) can be written as, 
\begin{equation}
\chi = \frac{\braket{{\mathbf{M}}^2}-\braket{{\mathbf{M}}}^2}{3 \epsilon_o \; k_B T V}
\end{equation}
\section{Boltzmann average of a quadrupole moment}
The average quadrupole moment at temperature T in the presence of uniform applied
field gradient is given by,
\begin{equation}
\braket{q_{\alpha\beta}} \;=\; \frac{\int d\Omega\; e^{\frac{-U(r)}{k_B T}}q_{\alpha\beta}}{\int d\Omega\; e^{\frac{-U(r)}{k_B T}}} \;=\; \frac{\int d\Omega\; e^{\frac{q_{\mu\nu}\;\partial_\nu E_\mu}{k_B T}}q_{\alpha\beta}}{\int d\Omega\; e^{\frac{q_{\mu\nu}\;\partial_\nu E_\mu}{k_B T}}},
\label{eq:1}
\end{equation}
where $\int \Omega = \int_0^{2\pi} \int_0^\pi \int_0^{2\pi}
sin\theta\; d\theta\ d\phi\ d\psi$ is the integration over Euler
angles and $ U(r) = -q_{\mu\nu}\;\partial_\nu E_\mu $ is the energy of
a quadrupole in the gradient of the  
applied field. The energy and quadrupole moment can be transformed
into body frame using following relation, 
\begin{equation}
\begin{split}
q_{\alpha\beta} &= \eta_{\alpha\alpha'}{q}^* _{\alpha'\beta'}\eta_{\beta'\beta} \\
U(r) &= - q:\vec{\nabla}\vec{E} = - q_{\mu\nu}\;\partial_\nu E_\mu = -\eta_{\mu\mu'}{q}^*_{\mu'\nu'}\eta_{\nu'\nu}\;\partial_\nu E_\mu.
\end{split}
\label{eq:2}
\end{equation}
Here the starred tensors are the components in the body fixed
frame. Substituting equation (\ref{eq:2}) in the equation (\ref{eq:1})
and taking linear terms in the expansion we get,
\begin{equation}
\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)},
\label{eq:3}
\end{equation}
where $\eta_{\alpha\alpha'}$ is the rotation matrix that transforms
the body fixed to the space fixed frame,
\[\eta_{\alpha\alpha'} 
= \left(\begin{array}{ccc}
cos\phi\; cos\psi - cos\theta\; sin\phi\; sin\psi & cos\theta\; cos\psi\; sin\phi + cos\phi\; sin\psi & sin\theta\; sin\phi \\
-cos\psi\; sin\phi - cos\theta\; cos\phi \; sin\psi & cos\theta\; cos\phi\; cos\psi - sin\phi\; sin\psi & cos\phi\; cos\theta \\
sin\theta\; sin\psi & -cos\psi\; sin\theta & cos\theta
\end{array} \right).\]
Integration of 1st and 2nd terms in the denominator gives $8 \pi^2$
and $8 \pi^2 /3\;\vec{\nabla}.\vec{E}\; Tr(q^*) $ respectively. The
second term vanishes for charge free space
(i.e. $\vec{\nabla}.\vec{E} \; = \; 0)$. Similarly integration of the
1st term in the numerator produces
$8 \pi^2 /3\; Tr(q^*)\delta_{\alpha\beta}$ and the 2nd term produces
$8 \pi^2 /15k_B T (3{q}^*_{\alpha'\beta'}{q}^*_{\beta'\alpha'} -
{q}^*_{\alpha'\alpha'}{q}^*_{\beta'\beta'})\partial_\alpha E_\beta$,
if $\vec{\nabla}.\vec{E} \; = \; 0$,
$ \partial_\alpha E_\beta = \partial_\beta E_\alpha$ and
${q}^*_{\alpha'\beta'}= {q}^*_{\beta'\alpha'}$. Therefore the
Boltzmann average of a quadrupole moment can be written as,

\begin{equation}
\braket{q_{\alpha\beta}}\; = \; \frac{1}{3} Tr(q^*)\;\delta_{\alpha\beta} + \frac{{\bar{q_o}}^2}{15k_BT}\;\partial_\alpha E_\beta,
\label{eq:4}
\end{equation}
here ${\bar{q_o}}^2=
3{q}^*_{\alpha'\beta'}{q}^*_{\beta'\alpha'}-{q}^*_{\alpha'\alpha'}{q}^*_{\beta'\beta'}$
which is a net quadrupole moment of a molecule.

To find the average quadrupole moment for the entire system, we can
multiply above equation with number of molecules $N$.
\begin{equation}
\braket{Q_{\alpha\beta}} \;=\; N\braket{q_{\alpha\beta}}\; = \; \frac{N}{3} Tr(q^*)\;\delta_{\alpha\beta} + \frac{{N\;\bar{q_o}}^2}{15k_BT}\;\partial_\alpha E_\beta 
\label{eq:5}
\end{equation}
The analogous relation in a system of the $N$ quadrupoles is,
\begin{equation}
\braket{Q_{\alpha\beta}} \; = \; \frac{1}{3}\;Tr(Q)\;\delta_{\alpha\beta} + \frac{{\braket{Q^2}}_{en} - {\braket{Q}}_{en}^2}{15k_BT}\;\partial_\alpha E_\beta 
\label{eq:6}
\end{equation}
where $\braket{{Q^2}}_{en}$ and $\braket{Q}_{en}^2$ are the ensemble
average of the square of the magnitude of the box quadrupole moment
and the square of the ensemble average of the box quadrupole moment
respectively. The ensemble averages in the above equation can be
obtained using relations,
${\braket{Q^2}}_{en} =
{\braket{3Q_{\alpha'\beta'}Q_{\beta'\alpha'}-Q_{\alpha'\alpha'}Q_{\beta'\beta'}}}_{en}$
and
${\braket{Q}}_{en}^2 = 3{\braket{Q_{\alpha'\beta'}}}_{en}\;
\braket{Q_{\beta'\alpha'}}_{en} -{\braket{Q_{\alpha'\alpha'}}}_{en}\;
{\braket{Q_{\beta'\beta'}}}_{en}$
respectively. The non-zero value of the $\braket{Q}_{en}^2$ is due to
the applied gradient, therefore its contribution is subtracted to find
the actual fluctuation in the box quadrupole moment.  From the
equations (\ref{eq:5}) and (\ref{eq:6}) we get,
\begin{equation}
\begin{split}
\braket{{Q^2}}_{en}-{\braket{Q}_{en}}^2 &= N\; {\bar{q_o}}^2\\
& and \\
\braket{Q_{\alpha\beta}}\; &= \; N\braket{q_{\alpha\beta}}.
\end{split}
\label{eq:7}
\end{equation}


\section{Relation between quadrupolar polarization and applied field gradient}
In the presence of the field gradient $\partial_\alpha {E}_\beta $, a
non-vanishing quadrupolar polarization (quadrupole moment per unit
volume) $\bar{Q}_{\alpha\beta}$ will be induced in the entire volume
of a sample. The total field at any point $\vec{r}$ in the sample is
given by,
\begin{equation}
\partial_\alpha E_\beta(\vec{r}) = \partial_\alpha {E^o}_\beta(\vec{r}) + \frac{1}{4\pi \epsilon_o}\int T_{\alpha\beta\gamma\delta}(\vec{r}-\vec{r'})\;\bar{Q}_{\gamma\delta}(\vec{r'})\; d^3r'
\label{eq:20}
\end{equation}
where $\partial_\alpha {E^o}_\beta$ is the applied field gradient and $ T_{\alpha\beta\gamma\delta}$ is the quadrupole-quadrupole interaction tensor which is obtained from applications of four successive gradient operator to the electrostatic kernal $v(r)$. 
\begin{equation}
\begin{split}
T_{\alpha\beta\gamma\delta}(r) &=\nabla_\alpha \nabla_\beta \nabla_\gamma \nabla_\delta\;v(r) \\ &= \left(\delta_{\alpha\beta}\delta_{\gamma\delta} + \delta_{\alpha\gamma}\delta_{\beta\delta}+ \delta_{\alpha\delta}\delta_{\beta\gamma}\right)v_{41}(r) + \left(\delta_{\alpha\beta} r_\gamma r_\delta + 5 \; permutations \right) \frac{v_{42}(r)}{r^2} \\ &+ r_\alpha r_\beta r_\gamma r_\delta\; \frac{v_{43}(r)}{r^4}, 
\end{split}
\label{quadRadial}
\end{equation}
where $v_{41}(r),\; v_{42}(r), \; and v_{43}(r)$ are defined in Paper I of the series (Refer Paper I) which have different functional forms for different cutoff methods. The integral in equation
(\ref{eq:20}) can be divided into two parts, $|\vec{r}-\vec{r'}|< R $
and $|\vec{r}-\vec{r'}|> R $ where $R \rightarrow 0 $. Since the total
field gradient vanishes at the singularity (see Appendix \ref{singular}),
equation (\ref{eq:20}) can be written as,
\begin{equation}
\partial_\alpha E_\beta(\vec{r}) = \partial_\alpha {E^o}_\beta(\vec{r}) +
  \frac{1}{4\pi \epsilon_o}\lim_{R \to
    0}\int\limits_{|\vec{r}-\vec{r'}|> R }
  T_{\alpha\beta\gamma\delta}(|\vec{r}-\vec{r'}|)\;\bar{Q}_{\gamma\delta}(\vec{r'})\;
  d^3r'. 
\label{eq:23}
\end{equation}
If $r = r'$ is excluded from the integration, the gradient of the electric can be expressed in terms of traceless quadrupole moment as below (Refer eq 2.7 of Logan I) ,
\begin{equation}
\partial_\alpha E_\beta(\vec{r}) = \partial_\alpha {E^o}_\beta(\vec{r}) + \frac{1}{12\pi \epsilon_o}\int\limits_{|\vec{r}-\vec{r'}|> R } T_{\alpha\beta\gamma\delta}(|\vec{r}-\vec{r'}|)\;\bar{\Theta}_{\gamma\delta}(\vec{r'})\; d^3r',
\label{eq:22}
\end{equation}
where $\Theta_{\alpha\beta} = 3Q_{\alpha\beta} - \delta_{\alpha\beta}Tr(Q)$
is the traceless quadrupole moment. The total quadrupolar polarization is written as,
\begin{equation}
\bar{Q}_{\alpha\beta}(\vec{r}) =  \frac{1}{3}\delta_{\alpha\beta}\;Tr(\bar{Q})+\epsilon_o {\chi}_Q\;\partial_\alpha E_\beta(\vec{r}),
\label{eq:21}
\end{equation}
where ${\chi}_Q = \frac{\braket{{Q^2}} - \braket{Q}^2}{15\epsilon_o Vk_BT}$ is the susceptibility of the quadrupolar fluid. In terms of traceless quadrupole moment, equation (\ref{eq:21}) can be written as,
\begin{equation}
\frac{1}{3}\bar{\Theta}_{\alpha\beta}(\vec{r}) =  \epsilon_o {\chi}_Q \; \partial_\alpha E_\beta (\vec{r})=  \epsilon_o {\chi}_Q \left(\partial_\alpha {E^o}_\beta(\vec{r}) + \frac{1}{12\pi \epsilon_o}\int\limits_{|\vec{r}-\vec{r'}|> R } T_{\alpha\beta\gamma\delta}(|\vec{r}-\vec{r'}|)\;\bar{\Theta}_{\gamma\delta}(\vec{r'})\; d^3r'\right)
\label{eq:22}
\end{equation}
For toroidal boundary conditions, both $\partial_\alpha E_\beta$ and
$\bar{\Theta}_{\alpha\beta}$ are uniform over the entire space. After
performing a Fourier transform (see the Appendix in the Neumann Paper),
we get,
\begin{equation}
\frac{1}{3}\widetilde{\bar{\Theta}}_{\alpha\beta}(\vec{k})=
\epsilon_o {\chi}_Q \;\left[{\partial_\alpha
    {E^o}_\beta}(\vec{k})+ \frac{1}{12\pi
    \epsilon_o}\;\widetilde{T}_{\alpha\beta\gamma\delta}(\vec{k})\;
 {\bar{\Theta}}_{\gamma\delta}(\vec{k})\right] 
\label{eq:24}
\end{equation} 
Since the quadrupolar polarization is in the direction of the applied
field
$\widetilde{\bar{\Theta}}_{\gamma\delta}(\vec{k}) =
\widetilde{\bar{\Theta}}_{\alpha\beta}(\vec{k})$
and
$\widetilde{T}_{\alpha\beta\gamma\delta}(\vec{k}) =
\widetilde{T}_{\alpha\beta\alpha\beta}(\vec{k})$.
(Note: This is true for orientational polarization but might not be
true for polarization due to changes in the atomic distribution. For
the orientational case, we need only to work with individual
components of the matrix.)
\begin{equation}
\begin{split}
\frac{1}{3}\widetilde{\bar{\Theta}}_{\alpha\beta}(\vec{k}) &= \epsilon_o \bar{\chi}_Q \left[{\partial_\alpha E^o_\beta}(\vec{k})+ \frac{1}{12\pi \epsilon_o}\widetilde{T}_{\alpha\beta\alpha\beta}(\vec{k})\;\widetilde{\bar{\Theta}}_{\alpha\beta}(\vec{k})\right] \\
&= \epsilon_o {\chi}_Q\;\left(1-\frac{1}{4\pi} {\chi}_Q
  \widetilde{T}_{\alpha\beta\alpha\beta}(\vec{k})\right)^{-1}
{\partial_\alpha E^o_\beta}(\vec{k}) 
\end{split}
\label{eq:25}
\end{equation} 
Matrix contraction of quadrupole-quadrupole tensor and quadrupole
moment gives the identity matrix. Therefore each component can be
treated as scalar. If the field gradient is homogeneous over the
entire volume, ${\partial_ \alpha E_\beta}(\vec{k}) = 0 $ except at
$ \vec{k} = 0$, hence it is sufficient to know
$ \widetilde{T}_{\alpha\beta\alpha\beta}(\vec{k})$ at $ \vec{k} =
0$. Therefore equation (\ref{eq:25}) can be written as,
\begin{equation}
\begin{split}
\frac{1}{3}\widetilde{\bar{\Theta}}_{\alpha\beta}({0}) &= \epsilon_o {\chi}_Q\left(1-\frac{1}{4\pi} {\chi}_Q\;\widetilde{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1} \partial_\alpha E^o_\beta({0})
\end{split}
\label{eq:26}
\end{equation}
where $ \widetilde{T}_{\alpha\beta\alpha\beta}(\vec{0})$ can be evaluated as,
\begin{equation}
\widetilde{T}_{\alpha\beta\alpha\beta}({0}) = \int \widetilde{T}_{\alpha\beta\alpha\beta}\;(\vec{r})d^3r
\label{eq:27}
\end{equation}
The susceptibility is modified as,
\begin{equation}
\Xi_Q = \bar{\chi}_Q\left(1-\frac{1}{4\pi} \bar{\chi}_Q\;\widetilde{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1}.
\label{eq:28}
\end{equation}
In terms of traced quadrupole moment equation (\ref{eq:25}) can be written as,
\begin{equation}
{\bar{Q}}_{\alpha\beta} = \frac{1}{3}\delta_{\alpha\beta}\;Tr(\bar{Q}) + \epsilon_o\; \Xi_Q\; \partial_\alpha E^o_\beta
\label{eq:29}
\end{equation}
From equations (\ref{eq:6}), (\ref{eq:28}), and (\ref{eq:29}) we get,
\begin{equation}
\begin{split}
&\frac{\braket{{Q^2}}_{en} - \braket{Q}_{en}^2}{15 \epsilon_o Vk_BT}\; =\; {\chi}_Q\left(1-\frac{1}{4\pi} {\chi}_Q\;\widetilde{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1}, \\ 
&{\chi}_Q \;=\; \frac{\braket{{Q^2}}_{en} - \braket{Q}_{en}^2}{15 \epsilon_o Vk_BT}\left(1 + \frac{1}{4\pi} \frac{\braket{{Q^2}}_{en} - \braket{Q}_{en}^2}{15 \epsilon_o Vk_BT}\;\widetilde{T}_{\alpha\beta\alpha\beta}({0})\right)^{-1}
\end{split}
\end{equation}
Finally in terms of quadrupolar correction factor ($A_{quad}$) the quadrupolar susceptibility is written as, 
\begin{equation}
{\chi}_Q \;=\; \frac{\braket{{Q^2}}_{en} - \braket{Q}_{en}^2}{15 \epsilon_o Vk_BT}\left(1 + \frac{\braket{{Q^2}}_{en} - \braket{Q}_{en}^2}{15 \epsilon_o Vk_BT}\; A_{quad}\right)^{-1}
\end{equation}
The $A_{quad}$ has dimension of the $length^{-2}$ is different for different cutoff methods which is listed in Table I. 
\section{The effective dielectric constant of a quadrupolar fluid}
The dielectric constant of a material is a measure of the screening of
the electrostatic energy when a material is placed in the
electrostatic field. If $U_{field}(r)$ and $U_{int}(r)$ are the total
electrostatic energies of a system due to the applied field and
interaction of the quadrupole with the gradient of the field
respectively. The dielectric constant can be written as,
\begin{equation}
\epsilon = \frac{U_{field} + U_{int}}{U_{field}}
\label{eq:31}
\end{equation}
The energy due to field can be written as,
\begin{equation}
U_{field}(r) = \int_V u(r)dv = \int_V \frac{1}{2}\epsilon_o {|E^o|}^2 dv,
\label{eq:32}
\end{equation}
where $u(r)$ is the energy density for applied electric field and $V$ is the total volume
of the system. The interaction energy can be assumed to be created by assembling quadrupoles from infinitely far away against the existing field gradient(Analogy Jackson's equations 4.83 - 4.94). Therefore
\begin{equation}
U_{int}(r) = -\frac{1}{2}\int_V \bar{Q}_{\alpha\beta}\;\partial_\beta E^o_{\alpha}\; dv,
\label{eq:33}
\end{equation}
where $\bar{Q}_{\alpha\beta}$ is the quadrupolar polarization. The quadrupoloarization can be expressed as linear function of the gradient of the applied field as below,
\begin{equation}
\bar{Q}_{\alpha\beta} = \epsilon_{o}\; {\chi_Q}\; \partial_\alpha E_\beta \approx \epsilon_{o}\; {\chi_Q}\; \partial_\alpha E^o_\beta ,
\label{eq:34}
\end{equation}
where ${\chi_Q} $ is a quadrupolar susceptibility. The local field gradient can be ignored in the region of low polarizability. In terms of the box quadrupole
moment it can be expressed as,
\begin{equation}
\chi_Q \;=\; \frac{\braket{{Q^2}}_{en} - \braket{Q}_{en}^2}{15 \epsilon_o Vk_BT}\left(1 + \frac{\braket{{Q^2}}_{en} - \braket{Q}_{en}^2}{15 \epsilon_o Vk_BT}\;A_{quad}\right
)^{-1}.
\label{eq:36}
\end{equation}
Using above relations (equation (\ref{eq:32}\;-\;\ref{eq:34}) in the
equation (\ref{eq:31}) produce,
\begin{equation}
\epsilon = 1 + {\chi}_Q \frac{\int_V |\partial_\alpha E^o_\beta|^2 dv}{\int_V{|E^o|}^2 dv}.
\end{equation}

\newpage

\appendix
\section{Point-multipolar interactions with a spatially-varying electric field}

We can treat objects $a$, $b$, and $c$ containing embedded collections
of charges. When we define the primitive moments, we sum over that
collections of charges using a local coordinate system within each
object.  The point charge, dipole, and quadrupole for object $a$ are
given by $C_a$, $\mathbf{D}_a$, and $\mathsf{Q}_a$, respectively.
These are the primitive multipoles which can be expressed as a
distribution of charges,
\begin{align}
C_a =&\sum_{k \, \text{in }a} q_k , \label{eq:charge} \\
D_{a\alpha} =&\sum_{k \, \text{in }a} q_k r_{k\alpha}, \label{eq:dipole}\\
Q_{a\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in }  a} q_k
r_{k\alpha}  r_{k\beta} . \label{eq:quadrupole}
\end{align}
Note that the definition of the primitive quadrupole here differs from
the standard traceless form, and contains an additional Taylor-series
based factor of $1/2$.  In Paper 1, we derived the forces and torques
each object exerts on the others.

Here we must also consider an external electric field that varies in
space: $\mathbf E(\mathbf r)$.  Each of the local charges $q_k$ in
object $a$ will then experience a slightly different field.  This
electric field can be expanded in a Taylor series around the local
origin of each object.  A different Taylor series expansion is carried
out for each object.

For a particular charge $q_k$, the electric field at that site's
position is given by:
\begin{equation}
E_\gamma + \nabla_\delta E_\gamma r_{k \delta} 
+ \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta}
r_{k \varepsilon} + ... 
\end{equation}
Note that the electric field is always evaluated at the origin of the
objects, and treating each object using point multipoles simplifies
this greatly.

To find the force exerted on object $a$ by the electric field, one
takes the electric field expression, and multiplies it by $q_k$, and
then sum over all charges in $a$:

\begin{align}
F_\gamma &=  \sum_{k \textrm{~in~} a} q_k \lbrace E_\gamma + \nabla_\delta E_\gamma r_{k \delta} 
+ \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma r_{k \delta}
r_{k \varepsilon} + ...  \rbrace  \\
 &= C_a E_\gamma + D_{a  \delta} \nabla_\delta E_\gamma 
+ Q_{a \delta \varepsilon} \nabla_\delta \nabla_\varepsilon E_\gamma +
... 
\end{align}

Similarly, the torque exerted by the field on $a$ can be expressed as
\begin{align}
\tau_\alpha &=  \sum_{k \textrm{~in~} a} (\mathbf r_k \times q_k \mathbf E)_\alpha \\
 & =  \sum_{k \textrm{~in~} a} \epsilon_{\alpha \beta \gamma} q_k
 r_{k\beta} E_\gamma(\mathbf r_k) \\
 & = \epsilon_{\alpha \beta \gamma} D_\beta E_\gamma 
+ 2 \epsilon_{\alpha \beta \gamma} Q_{\beta \delta} \nabla_\delta
E_\gamma + ...
\end{align}

The last term is essentially identical with form derived by Torres del
Castillo and M\'{e}ndez Garrido,\cite{Torres-del-Castillo:2006uo} although their derivation
utilized a traceless form of the quadrupole that is different than the
primitive definition in use here.  We note that the Levi-Civita symbol
can be eliminated by utilizing the matrix cross product in an
identical form as in Ref. \onlinecite{Smith98}:
\begin{equation}
\left[\mathsf{A} \times \mathsf{B}\right]_\alpha = \sum_\beta
\left[\mathsf{A}_{\alpha+1,\beta} \mathsf{B}_{\alpha+2,\beta}
  -\mathsf{A}_{\alpha+2,\beta} \mathsf{B}_{\alpha+1,\beta} 
\right]
\label{eq:matrixCross}
\end{equation}
where $\alpha+1$ and $\alpha+2$ are regarded as cyclic permuations of
the matrix indices.  In table \ref{tab:UFT} we give compact
expressions for how the multipole sites interact with an external
field that has exhibits spatial variations.

\begin{table}
\caption{Potential energy $(U)$, force $(\mathbf{F})$, and torque
  $(\mathbf{\tau})$ expressions for a multipolar site embedded in an
  electric field with spatial variations, $\mathbf{E}(\mathbf{r})$.
\label{tab:UFT}}
\begin{tabular}{r|ccc}
  & Charge & Dipole & Quadrupole \\ \hline
$U$ &  $C \phi(\mathbf{r})$ & $-\mathbf{D} \cdot \mathbf{E}(\mathbf{r})$ & $- \mathsf{Q}:\nabla \mathbf{E}(\mathbf{r})$ \\
$\mathbf{F}$ & $C \mathbf{E}(\mathbf{r})$ & $+\mathbf{D} \cdot \nabla \mathbf{E}(\mathbf{r})$ &  $+\mathsf{Q} : \nabla\nabla\mathbf{E}(\mathbf{r})$ \\
$\mathbf{\tau}$ & & $\mathbf{D} \times \mathbf{E}(\mathbf{r})$ & $+2 \mathsf{Q} \times \nabla \mathbf{E}(\mathbf{r})$
\end{tabular}
\end{table}
\section{Gradient of the field due to quadrupolar polarization}
\label{singular}
In this section, we will discuss the gradient of the field produced by
quadrupolar polarization. For this purpose, we consider a distribution
of charge ${\rho}(r)$ which gives rise to an electric field
$\vec{E}(r)$ and gradient of the field $\vec{\nabla} \vec{E}(r)$
throughout space. The total gradient of the electric field over volume
due to the all charges within the sphere of radius $R$ is given by
(cf. Jackson equation 4.14):
\begin{equation}
\int_{r<R} \vec{\nabla}\vec{E}\;d^3r = -\int_{r=R} R^2 \vec{E}\;\hat{n}\; d\Omega
\label{eq:8}
\end{equation}
where $d\Omega$ is the solid angle and $\hat{n}$ is the normal vector
of the surface of the sphere which is equal to
$sin[\theta]cos[\phi]\hat{x} + sin[\theta]sin[\phi]\hat{y} +
cos[\theta]\hat{z}$
in spherical coordinates.  For the charge density ${\rho}(r')$, the
total gradient of the electric field can be written as (cf. Jackson
equation 4.16),
\begin{equation}
\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
\label{eq:9}
\end{equation}
The radial function in the equation (\ref{eq:9}) can be expressed in
terms of spherical harmonics as (cf. Jackson equation 3.70),
\begin{equation}
\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)
\label{eq:10}
\end{equation}
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,
\begin{equation}
\begin{split}
\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 \\
 &= -\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
'
\end{split}
\label{eq:11}
\end{equation} 
The gradient of the product of radial function and spherical harmonics
is given by (cf. Arfken, p.811 eq. 16.94):
\begin{equation}
\begin{split}
\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
{\partial}{\partial r}+\frac{l}{r} \right]f(r)\; Y_{l, l-1, m}(\theta, \phi).
\end{split}
\label{eq:12}
\end{equation}
Using equation (\ref{eq:12}) we get,
\begin{equation}
\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}},
\label{eq:13}
\end{equation}
where $ Y_{l,l+1,m}(\theta, \phi)$ is the vector spherical harmonics
which can be expressed in terms of spherical harmonics as shown in
below (cf. Arfkan p.811),
\begin{equation}
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},
\label{eq:14}
\end{equation}
where $C(l+1,1,l|m_1,m_2,m)$ is a Clebsch-Gordan coefficient and
$\hat{e}_{m_2}$ is a spherical tensor of rank 1 which can be expressed
in terms of Cartesian coordinates,
\begin{equation}
{\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}}
\label{eq:15}
\end{equation} 
The normal vector $\hat{n} $ can be expressed in terms of spherical tensor of rank 1 as shown in below,
\begin{equation}
\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)
\label{eq:16}
\end{equation}
The surface integral of the product of $\hat{n}$ and
${Y_{l+1}}^{m_1}(\theta, \phi)$ gives,
\begin{equation}
\begin{split}
\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 \\
&=  \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 \\
&=   \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),
\end{split}
\label{eq:17}
\end{equation}
where ${Y_{l}}^{-m} = (-1)^m\;{{Y_{l}}^{m}}^* $ and
$ \int {{Y_{l}}^{m}}^*\;{Y_{l'}}^{m'}\;d\Omega =
\delta_{ll'}\delta_{mm'} $.
Non-vanishing values of equation \ref{eq:17} require $l = 0$,
therefore the value of $ m = 0 $. Since the values of $ m_1$ are -1,
1, and 0 then $m_2$ takes the values 1, -1, and 0, respectively
provided that $m = m_1 + m_2$.  Equation \ref{eq:11} can therefore be
modified,
\begin{equation}
\begin{split}
\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(
1, 1, 0|0,0,0)\;{\hat{e}_{0}}{\hat{e}_{0}} ]\; d^3r'.
\end{split}
\label{eq:18} 
\end{equation}
After substituting ${Y^*}_{00} = \frac{1}{\sqrt{4\pi}} $ and using the
values of the Clebsch-Gorden coefficients:  $ C(1, 1, 0|-1,1,0) =
\frac{1}{\sqrt{3}}, \;   C(1, 1, 0|-1,1,0)= \frac{1}{\sqrt{3}}$ and $
C(1, 1, 0|0,0,0) = -\frac{1}{\sqrt{3}}$ in equation \ref{eq:18} we
obtain,
\begin{equation}
\begin{split}
\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)\\
&= - \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).
\end{split}
\label{eq:19} 
\end{equation}
Equation (\ref{eq:19}) gives the total gradient of the field over a
sphere due to the distribution of the charges. For quadrupolar fluids
the total charge within a sphere is zero, therefore
$ \int_{r<R} \vec{\nabla}\vec{E}\;d^3r = 0 $.  Hence the quadrupolar
polarization produces zero net gradient of the field inside the
sphere.

\newpage

\bibliography{multipole}

\end{document}
Revision:	4353
Committed:	Tue Sep 15 15:56:07 2015 UTC (9 years, 10 months ago) by mlamichh
Content type:	application/x-tex
File size:	40716 byte(s)
Log Message:	Added new tex file for dielectric.