trunk/multipole/Dielectric_Supplemental.tex

% ****** Start of file aipsamp.tex ******
%
%   This file is part of the AIP files in the AIP distribution for REVTeX 4.
%   Version 4.1 of REVTeX, October 2009
%
%   Copyright (c) 2009 American Institute of Physics.
%
%   See the AIP README file for restrictions and more information.
%
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed
% as well as the rest of the prerequisites for REVTeX 4.1
%
% It also requires running BibTeX. The commands are as follows:
%
%  1)  latex  aipsamp
%  2)  bibtex aipsamp
%  3)  latex  aipsamp
%  4)  latex  aipsamp
%
% Use this file as a source of example code for your aip document.
% Use the file aiptemplate.tex as a template for your document.
\documentclass[%
 aip,jcp,
 amsmath,amssymb,
preprint,%
% reprint,%
%author-year,%
%author-numerical,%
jcp]{revtex4-1}

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


%\usepackage[mathlines]{lineno}% Enable numbering of text and display math
%\linenumbers\relax % Commence numbering lines

\begin{document}

\title{Supplemental Material for: 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

\maketitle

This document contains derivations of useful relationships for
electric field gradients and their interactions with point multipoles.
We rely heavily on both the notation and results from Torres del
Castillo and Mend\'{e}z Garido.\cite{Torres-del-Castillo:2006uo} In
this work, tensors are expressed in Cartesian components, using at
times a dyadic notation. This proves quite useful for our work as we
employ toroidal boundary conditions in our simulations, and these are
easily implemented in Cartesian coordinate systems.

An alternative formalism uses the theory of angular momentum and
spherical harmonics and is common in standard physics texts such as
Jackson,\cite{Jackson98} Morse and Feshbach, and Baym. Because this
approach has its own advantages, relationships are provided below
comparing that terminology to the Cartesian tensor notation.

The gradient of the electric field,
\begin{equation*}
\mathsf{G}(\mathbf{r}) = -\nabla \nabla \Phi(\mathbf{r}),
\end{equation*}
where $\Phi(\mathbf{r})$ is the electrostatic potential.  In a
charge-free region of space, $\nabla \cdot \mathbf{E}=0$, and
$\mathsf{G}$ is a symmetric traceless tensor.  From symmetry
arguments, we know that this tensor can be written in terms of just
five independent components.

Following Torres del Castillo and Mend\'{e}z Garido's notation, the
gradient of the electric field may also be written in terms of two
vectors $\mathbf{a}$ and $\mathbf{b}$,
\begin{equation*}
G_{ij}=\frac{1}{2} (a_i b_j + a_j b_i) - \frac{1}{3}(\mathbf a \cdot \mathbf b) \delta_{ij} .
\end{equation*} 
If the vectors $\mathbf{a}$ and $\mathbf{b}$ are unit vectors, the
electrostatic potential that generates a uniform gradient may be
written:
\begin{align}
\Phi(x, y, z) =\; -g_o & \left(\frac{1}{2}(a_1\;b_1 - \frac{cos\psi}{3})\;x^2+\frac{1}{2}(a_2\;b_2 - \frac{cos\psi}{3})\;y^2 + \frac{1}{2}(a_3\;b_3 - \frac{cos\psi}{3})\;z^2 \right. \\
 & \left. + \frac{(a_1\;b_2 + a_2\;b_1)}{2} x\;y + \frac{(a_1\;b_3 + a_3\;b_1)}{2} x\;z +  \frac{(a_2\;b_3 + a_3\;b_2)}{2} y\;z \right) .
\label{eq:appliedPotential}
\end{align} 
Note $\mathbf{a}\cdot\mathbf{a} = \mathbf{b} \cdot \mathbf{b} = 1$,
$\mathbf{a} \cdot \mathbf{b}=\cos \psi$, and $g_0$ is the overall
strength of the potential.  

An alternative to this notation is to write an electrostatic potential
that generates a uniform gradient using the notation of Morse and
Feshbach,
\begin{equation} \label{eq:quad_phi} 
\Phi(x,y,z) = - \left[ a_{20} \frac{2 z^2 -x^2 - y^2}{2}
+ 3 a_{21}^e \,xz + 3 a_{21}^o \,yz  
 + 6a_{22}^e \,xy +  3 a_{22}^o (x^2 - y^2) \right]  .
\end{equation}
Here we use the standard $(l,m)$ form for the $a_{lm}$ coefficients,
with superscript $e$ and $o$ denoting even and odd, respectively.
This form makes the functional analogy to ``d'' atomic states
apparent. The gradient of the electric field in this form is:
\begin{equation} \label{eq:grad_e2}
\mathsf{G} = 
\begin{pmatrix}
6 a_{22}^o - a_{20} & 6a_{22}^e & 3a_{21}^e\\
6a_{22}^e & -(a_{20}+6a_{22}^o) & 3a_{21}^o \\
3a_{21}^e  &  3a_{21}^o & 2a_{20} \\
\end{pmatrix} \\
\end{equation}
which can be factored as
\begin{gather}
\begin{aligned}
\mathsf{G}  = a_{20} 
\begin{pmatrix}
-1 & 0 & 0\\
0 & -1 & 0\\
0 & 0 & 2\\
\end{pmatrix}
+3a_{21}^e
\begin{pmatrix}
0 & 0 & 1\\
0 & 0 & 0\\
1 & 0 & 0\\
\end{pmatrix}
+3a_{21}^o
\begin{pmatrix}
0 & 0 & 0\\
0 & 0 & 1\\
0 & 1 & 0\\
\end{pmatrix}
+6a_{22}^e
\begin{pmatrix}
0 & 1 & 0\\
1 & 0 & 0\\
0 & 0 & 0\\
\end{pmatrix}
+6a_{22}^o
\begin{pmatrix}
1 & 0 & 0\\
0 & -1 & 0\\
0 & 0 & 0\\
\end{pmatrix}
\end{aligned}
 \label{eq:intro_tensors}
\end{gather}
The five matrices in the expression above represent five different
symmetric traceless tensors of rank 2.  The trace corresponds to
$\nabla \cdot \mathbf{E} = 0$, consistent with being in a charge-free
region.  Using the Cartesian notation of
Eq. (\ref{eq:appliedPotential}), this tensor is written:
\begin{equation}
\mathsf{G} =\nabla\bf{E} 
= \frac{g_o}{2}\left(\begin{array}{ccc}
2(a_1\; b_1 - \frac{cos\psi}{3}) &  (a_1\; b_2 \;+ a_2\; b_1) & (a_1\; b_3 \;+ a_3\; b_1) \\
 (a_2\; b_1 \;+ a_1\; b_2) & 2(a_2\; b_2 \;- \frac{cos\psi}{3}) & (a_2\; b_3 \;+ a_3\; b_3) \\
(a_3\; b_1 \;+ a_3\; b_2) & (a_3\; b_2 \;+ a_2\; b_3) & 2(a_3\; b_3 \;- \frac{cos\psi}{3})
\end{array} \right).
\label{eq:GC}
\end{equation}

It is useful to find vectors $\mathbf a$ and $\mathbf b$ that generate
the five types of tensors shown in Eq. (\ref{eq:intro_tensors}).  If
the two vectors are co-linear, e.g., $\psi=0$, $\mathbf{a}=(0,0,1)$ and
$\mathbf{b}=(0,0,1)$, then
\begin{equation*}
\mathsf{G} = \frac{g_0}{3}
\begin{pmatrix}
-1 & 0 & 0 \\
0 & -1 & 0 \\
0 & 0 & 2 \\
\end{pmatrix} ,
\end{equation*}
which is the $a_{20}$ symmetry.
To generate the $a_{22}^o$ symmetry, we take:
$\mathbf{a}= (\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}},0)$ and
$\mathbf{b}=(\frac{1}{\sqrt{2}}, -\frac{1}{\sqrt{2}},0)$
and find:
\begin{equation*}
\mathsf{G}=\frac{g_0}{2}
\begin{pmatrix}
1 & 0 & 0 \\
0 & -1 & 0 \\
0 & 0 & 0 \\
\end{pmatrix} .
\end{equation*}
To generate the $a_{22}^e$ symmetry, we take:
$\mathbf{a}= (1, 0, 0)$ and $\mathbf{b} = (0,1,0)$ and find:
\begin{equation*}
\mathsf{G}=\frac{g_0}{2}
\begin{pmatrix}
0 & 1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 0 \\
\end{pmatrix} .
\end{equation*}
The pattern is straightforward to continue for the other symmetries.

Using Eq. (\ref{eq:quad_phi}) the electric field is written:
\begin{equation}
\mathbf{E} 
= \left(\begin{array}{ccc}
\left(-a_{20} + 6a_{22}^o \right) x + 6a_{22}^e y + 3a_{21}^e  z  \\
6a_{22}^e x+(-a_{20} - 6a_{22}^o) y + 3a_{21}^e z \\
3a_{21}^e x +3a_{21}^o y + 2a_{20} z
\end{array} \right).
\label{eq:MFE}
\end{equation}
while using Eq. (\ref{eq:appliedPotential}), we find:
\begin{equation}
\mathbf{E}
=\frac{g_o}{2} \left(\begin{array}{ccc}
2(a_1\; b_1 - \frac{cos\psi}{3})\;x \;+  (a_1\; b_2 \;+ a_2\; b_1)\;y + (a_1\; b_3 \;+ a_3\; b_1)\;z \\
 (a_2\; b_1 \;+ a_1\; b_2)\;x + 2(a_2\; b_2 \;- \frac{cos\psi}{3})\;y +  (a_2\; b_3 \;+ a_3\; b_3)\;z \\
(a_3\; b_1 \;+ a_3\; b_2)\;x +  (a_3\; b_2 \;+ a_2\; b_3)\;y + 2(a_3\; b_3 \;- \frac{cos\psi}{3})\;z 
\end{array} \right).
\label{eq:CE}
\end{equation}
We find the notation of Ref. \onlinecite{Torres-del-Castillo:2006uo}
to be helpful when creating specific types of constant gradient
electric fields in simulations. For this reason,
Eqs. (\ref{eq:appliedPotential}), (\ref{eq:GC}), and (\ref{eq:CE}) are
used in our code.

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

This section develops formulas for the force and torque exerted by an
external electric field, $\mathbf{E}(\mathbf{r})$, on object
$a$. Object $a$ has an embedded collection of charges and in
simulations will represent a molecule, ion, or a coarse-grained
substructure. We describe the charge distributions using primitive
multipoles defined in Ref. \onlinecite{PaperI} by
\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}
where $\mathbf{r}_k$ is the local coordinate system for the object
(usually the center of mass of object $a$).  Components of vectors and
tensors are given using Green indices, using the Einstein repeated
summation notation. 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 Ref.
\onlinecite{PaperI}, we derived the forces and torques each object
exerts on the other objects in the system.

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. For a particular charge $q_k$, the electric
field at that site's position is given by:
\begin{equation}
\mathbf{E}(\mathbf{r}_k) = E_\gamma|_{\mathbf{r}_k = 0} + \nabla_\delta E_\gamma |_{\mathbf{r}_k = 0}  r_{k \delta} 
+ \frac {1}{2} \nabla_\delta \nabla_\varepsilon E_\gamma|_{\mathbf{r}_k = 0}  r_{k \delta}
r_{k \varepsilon} + ... 
\end{equation}
Note that once one shrinks object $a$ to point size, the ${E}_\gamma$
terms are all evaluated at the center of the object (now a
point). Thus later the ${E}_\gamma$ terms can be written using the
same global origin for all objects $a, b, c, ...$ in the system. The
force exerted on object $a$ by the electric field is given by,

\begin{align}
F^a_\gamma = \sum_{k \textrm{~in~} a} E_\gamma(\mathbf{r}_k) &=  \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}
Thus in terms of the global origin $\mathbf{r}$, ${F}_\gamma(\mathbf{r}) = C {E}_\gamma(\mathbf{r})$ etc.
  
Similarly, the torque exerted by the field on $a$ can be expressed as
\begin{align}
\tau^a_\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}
We note that the Levi-Civita symbol can be eliminated by utilizing the matrix cross product as defined 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. Finally, the interaction energy $U^a$ of object $a$ with the external field is given by,
\begin{equation}
U^a = \sum_{k~in~a} q_k \phi_k (\mathrm{r}_k)
\end{equation}
Performing another Taylor series expansion about the local body origin,
\begin{equation}
\phi({\mathbf{r}_k}) = \phi|_{\mathbf{r}_k = 0 } + r_{k \alpha} \nabla_\alpha \phi_\alpha|_{\mathbf{r}_k = 0 } + \frac{1}{2} r_{k\alpha}r_{k\beta}\nabla_\alpha \nabla_\beta \phi|_{\mathbf{r}_k = 0} + ...
\end{equation}
Writing this in terms of the global origin $\mathrm{r}$, we find
\begin{equation}
U(\mathbf{r}) = \mathrm{C} \phi(\mathbf{r}) - \mathrm{D}_\alpha \mathrm{E}_\alpha - \mathrm{Q}_{\alpha\beta}\nabla_\alpha \mathrm{E}_\beta + ...
\end{equation}
The results has been summarized in Table I.

\begin{table}
\caption{Potential energy $(U)$, force $(\mathbf{F})$, and torque
  $(\mathbf{\tau})$ expressions for a multipolar site $\mathrm{r}$ in an
  electric field, $\mathbf{E}(\mathbf{r})$.
\label{tab:UFT}}
\begin{tabular}{r|ccc}
  & Charge & Dipole & Quadrupole \\ \hline
$U(\mathbf{r})$ &  $C \phi(\mathbf{r})$ & $-\mathbf{D} \cdot \mathbf{E}(\mathbf{r})$ & $- \mathsf{Q}:\nabla \mathbf{E}(\mathbf{r})$ \\
$\mathbf{F}(\mathbf{r})$ & $C \mathbf{E}(\mathbf{r})$ & $\mathbf{D} \cdot \nabla \mathbf{E}(\mathbf{r})$ &  $\mathsf{Q} : \nabla\nabla\mathbf{E}(\mathbf{r})$ \\
$\mathbf{\tau}(\mathbf{r})$ & & $\mathbf{D} \times \mathbf{E}(\mathbf{r})$ & $2 \mathsf{Q} \times \nabla \mathbf{E}(\mathbf{r})$
\end{tabular}
\end{table}


\section{Boltzmann averages for orientational polarization}
The dielectric properties of the system is mainly arise from two
different ways: i) the applied field distort the charge distributions
so it produces an induced multipolar moment in each molecule; and ii)
the applied field tends to line up originally randomly oriented
molecular moment towards the direction of the applied field. In this
study, we basically focus on the orientational contribution in the
dielectric properties. If we consider a system of molecules in the
presence of external field perturbation, the perturbation experienced
by any molecule will not be only due to external field or field
gradient but also due to the field or field gradient produced by the
all other molecules in the system. In the following subsections
\ref{subsec:boltzAverage-Dipole} and \ref{subsec:boltzAverage-Quad},
we will discuss about the molecular polarization only due to external
field perturbation. The contribution of the field or field gradient
due to all other molecules will be taken into account while
calculating correction factor in the paper.

\subsection{Dipoles}
\label{subsec:boltzAverage-Dipole}
Consider a system of molecules, each with permanent dipole moment
$p_o$. In the absense of external field, thermal agitation orients the
dipoles randomly, reducing the system moment to zero.  External fields
will tend to line up the dipoles in the direction of applied field.
Here we have considered net field from all other molecules is
considered to be zero.  Therefore the total Hamiltonian of each
molecule is,\cite{Jackson98}
\begin{equation}
H = H_o - \mathbf{p_o}\cdot \mathbf{E},
\end{equation}
where $H_o$ is a function of the internal coordinates of the molecule.
The Boltzmann average of the dipole moment is given by,
\begin{equation}
\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}}},
\end{equation}
where $\bf{E}$ is selected along z-axis. If we consider that the
applied field is small, \textit{i.e.} $\frac{p_oE\; cos\theta}{k_B T} << 1$, 
\begin{equation}
\braket{p_{mol}}  \approx \frac{1}{3}\frac{{p_o}^2}{k_B T}E,
\end{equation}
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.

\subsection{Quadrupoles}
\label{subsec:boltzAverage-Quad}
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}
\begin{equation}
\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}}},
\label{boltzQuad}
\end{equation}
where $\int d\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, $ H = H_o -q_{\mu\nu}\;\partial_\nu E_\mu $ is the energy of
a quadrupole in the gradient of the  
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, 
\begin{equation}
\begin{split}
&q_{\alpha\beta} = \eta_{\alpha\alpha'}\;\eta_{\beta\beta'}\;{q}^* _{\alpha'\beta'} \\
&H = H_o - q:{\nabla}\mathbf{E} = H_o - q_{\mu\nu}\;\partial_\nu E_\mu = H_o -\eta_{\mu\mu'}\;\eta_{\nu\nu'}\;{q}^*_{\mu'\nu'}\;\partial_\nu E_\mu.
\end{split}
\label{energyQuad}
\end{equation}
Here the starred tensors are the components in the body fixed
frame. Substituting equation (\ref{energyQuad}) in the equation (\ref{boltzQuad})
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'}\;\eta_{\nu\nu'}\;{q}^*_{\mu'\nu'}\;\partial_\nu E_\mu }{k_B T}\right)q_{\alpha\beta}}{ \int d\Omega \left(1 + \frac{\eta_{\mu\mu'}\;\eta_{\nu\nu'}\;{q}^*_{\mu'\nu'}\;\partial_\nu E_\mu }{k_B T}\right)},
\end{equation}
where $\eta_{\alpha\alpha'}$ is the inverse of the rotation matrix that transforms
the body fixed co-ordinates to the space co-ordinates,
\[\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\; sin\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\;{\nabla}.\mathbf{E}\; Tr(q^*) $ respectively. The
second term vanishes for charge free space, ${\nabla}.\mathbf{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 ${\nabla}.\mathbf{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,
\end{equation}
where $ \alpha_q = \frac{{\bar{q_o}}^2}{15k_BT} $ is a molecular quadrupole polarizablity  and  ${\bar{q_o}}^2=
3{q}^*_{\alpha'\beta'}{q}^*_{\beta'\alpha'}-{q}^*_{\alpha'\alpha'}{q}^*_{\beta'\beta'}$ is a square of the net quadrupole moment of a molecule. 

% \section{External application of a uniform field gradient}
% \label{Ap:fieldOrGradient}

% To satisfy the condition $ \nabla \cdot \mathbf{E} = 0 $, within the box of molecules we have taken electrostatic potential in the following form
% \begin{equation}
% \begin{split}
% \phi(x, y, z) =\; &-g_o \left(\frac{1}{2}(a_1\;b_1 - \frac{cos\psi}{3})\;x^2+\frac{1}{2}(a_2\;b_2 - \frac{cos\psi}{3})\;y^2 + \frac{1}{2}(a_3\;b_3 - \frac{cos\psi}{3})\;z^2 \right. \\
% & \left. + \frac{(a_1\;b_2 + a_2\;b_1)}{2} x\;y + \frac{(a_1\;b_3 + a_3\;b_1)}{2} x\;z +  \frac{(a_2\;b_3 + a_3\;b_2)}{2} y\;z \right),
% \end{split}
% \label{eq:appliedPotential}
% \end{equation} 
% where $a = (a_1, a_2, a_3)$ and $b = (b_1, b_2, b_3)$ are basis vectors  determine coefficients in x, y, and z direction. And $g_o$ and $\psi$ are overall strength of the potential and angle between basis vectors respectively. The electric field derived from the above potential is,
% \[\mathbf{E} 
% = \frac{g_o}{2} \left(\begin{array}{ccc}
% 2(a_1\; b_1 - \frac{cos\psi}{3})\;x \;+  (a_1\; b_2 \;+ a_2\; b_1)\;y + (a_1\; b_3 \;+ a_3\; b_1)\;z \\
%  (a_2\; b_1 \;+ a_1\; b_2)\;x + 2(a_2\; b_2 \;- \frac{cos\psi}{3})\;y +  (a_2\; b_3 \;+ a_3\; b_2)\;z \\
% (a_3\; b_1 \;+ a_3\; b_2)\;x +  (a_3\; b_2 \;+ a_2\; b_3)y + 2(a_3\; b_3 \;- \frac{cos\psi}{3})\;z 
% \end{array} \right).\] 
% The gradient of the applied field derived from the potential can be written in the following form,
% \[\nabla\mathbf{E} 
% = \frac{g_o}{2}\left(\begin{array}{ccc}
% 2(a_1\; b_1 - \frac{cos\psi}{3}) &  (a_1\; b_2 \;+ a_2\; b_1) & (a_1\; b_3 \;+ a_3\; b_1) \\
%  (a_2\; b_1 \;+ a_1\; b_2) & 2(a_2\; b_2 \;- \frac{cos\psi}{3}) & (a_2\; b_3 \;+ a_3\; b_2) \\
% (a_3\; b_1 \;+ a_3\; b_2) & (a_3\; b_2 \;+ a_2\; b_3) & 2(a_3\; b_3 \;- \frac{cos\psi}{3})
% \end{array} \right).\] 


\section{Gradient of the field due to quadrupolar polarization}
\label{singularQuad}
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}(\mathbf r)$ which gives rise to an electric field
$\mathbf{E}(\mathbf r)$ and gradient of the field ${\nabla} \mathbf{E}(\mathbf 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} {\nabla}\mathbf{E}\;d^3r = -\int_{r=R} R^2 \mathbf{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, 
$\hat{n} = sin[\theta]cos[\phi]\hat{x} + sin[\theta]sin[\phi]\hat{y} +
cos[\theta]\hat{z}$
in spherical coordinates.  For the charge density ${\rho}(\mathbf r')$, the
total gradient of the electric field can be written as, ~\cite{Jackson98}
\begin{equation}
\int_{r<R} {\nabla}\mathbf {E}\; d^3r=-\int_{r=R} R^2\; {\nabla}\Phi\; \hat{n}\; d\Omega  =-\frac{1}{4\pi\;\epsilon_o}\int_{r=R} R^2\; {\nabla}\;\left(\int \frac{\rho(\mathbf r')}{|\mathbf{r}-\mathbf{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,\cite{Jackson98}
\begin{equation}
\frac{1}{|\mathbf{r} - \mathbf{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} {\nabla}\mathbf{E}\;d^3r &=-\frac{R^2}{\epsilon_o}\int_{r=R} \; {\nabla}\;\left(\int \rho(\mathbf 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(\mathbf 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:\cite{Arfkan}
\begin{equation}
\begin{split}
{\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}
{\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 a vector spherical harmonics \cite{Arfkan}. Using Clebsch-Gorden coefficients $C(l+1, 1, l|m_1,m_2,m) $, equation \ref{eq:14} can be written in spherical harmonics,
\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}
Here $\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} $ is then 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} {\nabla}\mathbf{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' \\
&= -\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}
In the last step, the charge density was integrated over the sphere, yielding a total charge $\mathrm{C_total}$.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} {\nabla}\mathbf{E}\;d^3r = 0 $.  Hence the quadrupolar
polarization produces zero net gradient of the field inside the
sphere.


\bibliography{dielectric_new}
\end{document}
%
% ****** End of file multipole.tex ******
Revision:	4417
Committed:	Sun Apr 10 12:57:41 2016 UTC (8 years, 2 months ago) by gezelter
Content type:	application/x-tex
File size:	27648 byte(s)
Log Message:	Latest changes