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{array}
\newcolumntype{L}[1]{>{\raggedright\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
\newcolumntype{C}[1]{>{\centering\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
\newcolumntype{R}[1]{>{\raggedleft\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}

\renewcommand{\theequation}{S\arabic{equation}}
\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\thetable}{S\arabic{table}}
\renewcommand{\bibnumfmt}[1]{[S#1]~}
\renewcommand{\citenumfont}[1]{S#1}

%\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

\begin{abstract}
  This document includes useful relationships for computing the
  interactions between fields and field gradients and point multipolar
  representations of molecular electrostatics. We also provide
  explanatory derivations of a number of relationships used in the
  main text. This includes the Boltzmann averages of quadrupole
  orientations, and the interaction of a quadrupole density with the
  self-generated field gradient. This last relationship is assumed to
  be zero in the main text but is explicitly shown to be zero here.  A
  discussion of method-dependent corrections to the distance-dependent
  Kirkwood factors is also included.
\end{abstract}

\maketitle

\section{Generating Uniform Field Gradients}
One important task in carrying out the simulations mentioned in the
main text was to generate uniform electric field gradients.  To do
this, we relied 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 were expressed in Cartesian components, using at
times a dyadic notation. This proves quite useful for computer
simulations that make use of toroidal boundary conditions.

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,\cite{Morse:1946zr} and
Stone.\cite{Stone:1997ly} 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) =\; -\frac{g_o}{2} & \left(\left(a_1b_1 -
                         \frac{cos\psi}{3}\right)\;x^2+\left(a_2b_2
                         - \frac{cos\psi}{3}\right)\;y^2 +
                         \left(a_3b_3 -
                         \frac{cos\psi}{3}\right)\;z^2 \right. \nonumber \\
 & + (a_1b_2 + a_2b_1)\; xy + (a_1b_3 + a_3b_1)\; xz + (a_2b_3 + a_3b_2)\; yz \bigg) .
\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.

Taking the gradient of Eq. (\ref{eq:appliedPotential}), we find the
field due to this potential,
\begin{equation}
\mathbf{E} = -\nabla \Phi
=\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} 
while the gradient of the electric field in this form,
\begin{equation}
\mathsf{G} = \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_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}  
is uniform over the entire space.  Therefore, to describe a uniform
gradient in this notation, two unit vectors ($\mathbf{a}$ and
$\mathbf{b}$) as well as a potential strength, $g_0$, must be
specified. As expected, this requires five independent parameters.

The common alternative to the Cartesian notation expresses the
electrostatic potential using the notation of Morse and
Feshbach,\cite{Morse:1946zr}
\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. 

Applying the gradient operator to Eq. (\ref{eq:quad_phi}) the electric
field due to this potential,
\begin{equation}
\mathbf{E} = -\nabla \Phi
= \left(\begin{array}{ccc}
\left( 6a_{22}^o -a_{20} \right)\; x &+\; 6a_{22}^e\; y &+\; 3a_{21}^e\;  z  \\
6a_{22}^e\; x & -\; (a_{20} + 6a_{22}^o)\; y & +\; 3a_{21}^o\; z \\
3a_{21}^e\; x & +\; 3a_{21}^o\; y & +\; 2a_{20}\; z
\end{array} \right),
\label{eq:MFE}
\end{equation}
while 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 is also uniform over the entire space.  This form for the
gradient 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. 

It is useful to find the Cartesian 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.

We find the notation of Ref. \onlinecite{Torres-del-Castillo:2006uo}
helpful when creating specific types of constant gradient electric
fields in simulations. For this reason,
Eqs. (\ref{eq:appliedPotential}), (\ref{eq:CE}), and (\ref{eq:GC}) are
implemented in our code.  In the simulations using constant applied
gradients that are mentioned in the main text, we utilized a field
with the $a_{22}^e$ symmetry using vectors, $\mathbf{a}= (1, 0, 0)$
and $\mathbf{b} = (0,1,0)$.

\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$.\cite{Raab:2004ve} 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 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 if one shrinks object $a$ to a single point, 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 (molecular) origin for all point charges in the object. The force
exerted on object $a$ by the electric field is given by,
\begin{align}
F^a_\gamma = \sum_{k \textrm{~in~} a} q_k 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 $\mathbf{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}
These results have been summarized in Table \ref{tab:UFT}.

\begin{table}
\caption{Potential energy $(U)$, force $(\mathbf{F})$, and torque
  $(\mathbf{\tau})$ expressions for a multipolar site at $\mathbf{r}$ in an
  electric field, $\mathbf{E}(\mathbf{r})$ using the definitions of the multipoles in Eqs. (\ref{eq:charge}), (\ref{eq:dipole}) and (\ref{eq:quadrupole}).  
  \label{tab:UFT}}
\begin{tabular}{r|C{3cm}C{3cm}C{3cm}}
  & 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}
If we consider a collection of molecules in the presence of external
field, the perturbation experienced by any one molecule will include
contributions to the field or field gradient produced by the all other
molecules in the system. In subsections
\ref{subsec:boltzAverage-Dipole} and \ref{subsec:boltzAverage-Quad},
we discuss the molecular polarization due solely to external field
perturbations.  This illustrates the origins of the polarizability
equations (Eqs. 6, 20, and 21) in the main text.

\subsection{Dipoles}
\label{subsec:boltzAverage-Dipole}
Consider a system of molecules, each with permanent dipole moment
$p_o$. In the absense of an external field, thermal agitation orients
the dipoles randomly, and the system moment, $\mathbf{P}$, is zero.
External fields will line up the dipoles in the direction of applied
field.  Here we consider the net field from all other molecules to be
zero.  Therefore the total Hamiltonian acting on 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 in the direction of the
field is given by,
\begin{equation}
\langle p_{mol} \rangle = \frac{\displaystyle\int p_o \cos\theta
  e^{~p_o E \cos\theta /k_B T}\; d\Omega}{\displaystyle\int  e^{~p_o E \cos\theta/k_B
    T}\; d\Omega},
\end{equation}
where the $z$-axis is taken in the direction of the applied field,
$\bf{E}$ and
$\int d\Omega = \int_0^\pi \sin\theta\; d\theta \int_0^{2\pi} d\phi
\int_0^{2\pi} d\psi$
is an integration over Euler angles describing the orientation of the
molecule.

If the external fields are small, \textit{i.e.}
$p_oE \cos\theta / k_B T << 1$,
\begin{equation}
\langle p_{mol} \rangle \approx \frac{{p_o}^2}{3 k_B T}E,
\end{equation}
where $ \alpha_p = \frac{{p_o}^2}{3 k_B T}$ is the molecular
polarizability. The orientational polarization depends inversely on
the temperature as the applied field must overcome thermal agitation
to orient the dipoles.

\subsection{Quadrupoles}
\label{subsec:boltzAverage-Quad}
If instead, our system consists of molecules with permanent
\textit{quadrupole} tensor $q_{\alpha\beta}$. The average quadrupole
at temperature $T$ in the presence of uniform applied field gradient
is given by,\cite{AduGyamfi78, AduGyamfi81}
\begin{equation}
\langle q_{\alpha\beta} \rangle \;=\; \frac{\displaystyle\int
  q_{\alpha\beta}\; e^{-H/k_B T}\; d\Omega}{\displaystyle\int
  e^{-H/k_B T}\; d\Omega} \;=\; \frac{\displaystyle\int
  q_{\alpha\beta}\; e^{~q_{\mu\nu}\;\partial_\nu E_\mu /k_B T}\;
  d\Omega}{\displaystyle\int  e^{~q_{\mu\nu}\;\partial_\nu E_\mu /k_B
    T}\; d\Omega },
\label{boltzQuad}
\end{equation}
where $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 the body frame using a
rotation matrix $\mathsf{\eta}^{-1}$,
\begin{align}
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. \label{energyQuad}
\end{align}
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 obtain,
\begin{equation}
\braket{q_{\alpha\beta}} = \frac{\displaystyle \int q_{\alpha\beta} \left(1 +
    \frac{\eta_{\mu\mu'}\;\eta_{\nu\nu'}\;{q}^*_{\mu'\nu'}\;\partial_\nu
      E_\mu }{k_B T}\right)\;  d\Omega}{\displaystyle \int \left(1 + \frac{\eta_{\mu\mu'}\;\eta_{\nu\nu'}\;{q}^*_{\mu'\nu'}\;\partial_\nu E_\mu }{k_B T}\right)\; d\Omega}.
\end{equation}
Recall that $\eta_{\alpha\alpha'}$ is the inverse of the rotation
matrix that transforms the body fixed coordinates to the space
coordinates.
% \[\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 the first and second terms in the denominator gives
$8 \pi^2$ and
$8 \pi^2 ({\nabla} \cdot \mathbf{E}) \mathrm{Tr}(q^*) / 3 $
respectively. The second term vanishes for charge free space (where
${\nabla} \cdot \mathbf{E}=0$). Similarly, integration of the first
term in the numerator produces
$8 \pi^2 \delta_{\alpha\beta} \mathrm{Tr}(q^*) / 3$ while the second
produces
$8 \pi^2 (3{q}^*_{\alpha'\beta'}{q}^*_{\beta'\alpha'} -
{q}^*_{\alpha'\alpha'}{q}^*_{\beta'\beta'})\partial_\alpha E_\beta /
15 k_B T $.
Therefore the Boltzmann average of a quadrupole moment can be written
as,
\begin{equation}
\langle q_{\alpha\beta} \rangle =  \frac{1}{3} \mathrm{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 the square of the net quadrupole moment of a molecule.

\section{Gradient of the field due to quadrupolar polarization}
\label{singularQuad}
In section IV.C of the main text, we stated that for quadrupolar
fluids, the self-contribution to the field gradient vanishes at the
singularity. In this section, we prove this statement.  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
gradient of the electric field over volume due to the charges within
the sphere of radius $R$ is given by (cf. Ref. \onlinecite{Jackson98},
equation 4.14):
\begin{equation}
\int_{r<R} \nabla\mathbf{E} d\mathbf{r} = -\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, 
\begin{equation}
\hat{n} = \sin\theta\cos\phi\; \hat{x} + \sin\theta\sin\phi\; \hat{y} +
\cos\theta\; \hat{z}
\end{equation}
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\mathbf{r}=-\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\mathbf{r}'\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\mathbf{r} &=-\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\mathbf{r}'\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\mathbf{r}
' .
\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}
where $Y_{l,l+1,m}(\theta, \phi)$ is a vector spherical
harmonic.\cite{Arfkan} 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}
Using Clebsch-Gordan coefficients $C(l+1,1,l|m_1,m_2,m)$, the vector
spherical harmonics can be written in terms of 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\mathbf{r} = &- \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\mathbf{r}' \\
&= -\sqrt{\frac{4\pi}{{3}}}\;\frac{1}{\epsilon_o}\int \rho(r')\;d\mathbf{r}'\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_\mathrm{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 $C_\mathrm{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\mathbf{r} = 0 $.  Hence the quadrupolar
polarization produces zero net gradient of the field inside the
sphere.

\section{Corrections to the distance-dependent Kirkwood function}
\label{kirkwood}
In the main text, we provide data on the distance-dependent Kirkwood
function, 
\begin{equation}
G_K(r) = \left< \frac{1}{N} \sum_{i} \sum_{\substack{j \\ r_{ij} < r}}
  \frac{\mathbf{D}_i \cdot \mathbf{D}_j}{\left| D_i \right| \left| D_j
    \right|} \right>,
\label{eq:kirkwood}
\end{equation}
which is a sensitive measure of orientational ordering in dipolar
liquids.  We noted in the main text that because the Kirkwood function
is measuring the same bulk dipolar response as the dielectric
constant, it is possible to apply a correction for the truncation,
force shifting, and damping to arrive at a corrected Kirkwood
function. Starting with Eq. (16) in the main text and recognizing
$\epsilon = 1 + \chi_D$ and $\epsilon_{CB} = 1 + \alpha_D$, we arrive
at a relationship between the (corrected) susceptibility,
\begin{equation}
\chi_D = \frac{\alpha_D}{1 + (A-1) \frac{\alpha_D}{3}}.
\label{eq:chiAlpha}
\end{equation}
and the dipole polarizability, which is easily obtained from bulk
simulations, 
\begin{equation}
\alpha_D =\frac{\braket{\mathbf{M}^2}-{\braket{\mathbf{M}}}^2}{3
  \epsilon_o V k_B T}.
\end{equation}
In the absence of bulk dipolar ordering, ${\braket{\mathbf{M}}}^2=0$,
and we may recognize the polarizability as a limit of the
distance-dependent Kirkwood function,
\begin{equation}
\alpha_D = \left( \frac{\rho D^2}{3 \epsilon_0 k_B T}\right) \lim_{r->\infty}
G_K(r).
\label{eq:alphaGK}
\end{equation}
Here we have used the number density, $\rho = N/V$ and the molecular
dipole moment, $D$, to connect with Eq. (\ref{eq:kirkwood}).  By
analogy, the dipolar susceptibility may be connected with a corrected
(but unknown) version of the Kirkwood function,
\begin{equation}
\chi_D = \left(\frac{\rho D^2}{3 \epsilon_0 k_B T}\right) \lim_{r->\infty}
G_K^c(r).
\label{eq:chiGKc}
\end{equation}
Substituting the corrected and raw Kirkwood expressions,
Eqs. (\ref{eq:chiGKc}) and (\ref{eq:alphaGK}), into
Eq. (\ref{eq:chiAlpha}), and removing the limits, we obtain a
correction formula for the Kirkwood function,
\begin{equation}
G_K^c(r) = \frac{G_K(r)}{ 1 + (A-1) \frac{\rho D^2 G_K(r)}{9
    \epsilon_0 k_B T}}.
\label{eq:kirkwoodCorr}
\end{equation}
Note that this correction forumla uses the same $A$ parameter from
Table I in the main text.  As in case of the dielectric constant, the
Kirkwood correction is similarly sensitive to values of $A$ away from
unity.  In Fig.~\ref{fig:kirkwoodCorr} we show the corrected
$G_K^c(r)$ functions for the three real space methods with
$r_c = 3.52 \sigma = 12$~\AA~ and for the Ewald sum (with
$\kappa = 0.3119$~\AA$^{-1}$).

\begin{figure}
\includegraphics[width=\linewidth]{kirkwoodCorr.eps}
\caption{Corrected distance-dependent factors of the dipolar system
  for the three real space methods at a range of Gaussian damping
  parameters ($\alpha$) with a cutoff
  $r_c = 3.52 \sigma = 12$~\AA$^{-1}$.}
\label{fig:kirkwoodCorr}
\end{figure}

The correction does help reduce the effect of the ``hole'' in the
underdamped cases, particularly for the GSF method.  However, this
comes at the expense of a divergence when $G_K(r) \sim 1$ for the
underdamped SP case.  We find it more useful to look at the
uncorrected $G_K(r)$ functions to study orientational correlations,
which is why the uncorrected functions are shown in the main text.

\bibliography{dielectric_new}
\end{document}
%
% ****** End of file multipole.tex ******
Revision:	4434
Committed:	Mon Aug 1 13:21:53 2016 UTC (9 years, 1 month ago) by gezelter
Content type:	application/x-tex
File size:	30898 byte(s)
Log Message:	Latest changes for review 2