trunk/ripple2/ripple.tex

\documentclass[aps,pre,endfloats*,preprint,amssymb,showpacs]{revtex4}
\usepackage{epsfig}

\begin{document}
\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}

%\bibliographystyle{aps}

\title{Spontaneous Corrugation of Dipolar Membranes}
\author{Xiuquan Sun and J. Daniel Gezelter}
\email[]{E-mail: gezelter@nd.edu}
\affiliation{Department of Chemistry and Biochemistry,\\
University of Notre Dame, \\
Notre Dame, Indiana 46556}

\date{\today} 

\begin{abstract}
We present a simple model for dipolar membranes that gives
lattice-bound point dipoles complete orientational freedom as well as
translational freedom along one coordinate (out of the plane of the
membrane).  There is an additional harmonic surface tension which
binds each of the dipoles to the six nearest neighbors on either
hexagonal or distorted-hexagonal lattices.  The translational freedom
of the dipoles allows hexagonal lattices to find states that break out
of the normal orientational disorder of frustrated configurations and
which are stabilized by long-range antiferroelectric ordering.  In
order to break out of the frustrated states, the dipolar membranes
form corrugated or ``rippled'' phases that make the lattices
effectively non-hexagonal.  We observe three common features of the
corrugated dipolar membranes: 1) the corrugated phases develop easily
when hosted on hexagonal lattices, 2) the wave vectors for the surface
ripples are always found to be perpendicular to the dipole director
axis, and 3) on hexagonal lattices, the dipole director axis is found
to be parallel to any of the three equivalent lattice directions.
\end{abstract}

\pacs{68.03.Hj, 82.20.Wt}
\maketitle


\section{Introduction}
\label{Int}
There has been intense recent interest in the phase behavior of
dipolar
fluids.\cite{Tlusty00,Teixeira00,Tavares02,Duncan04,Holm05,Duncan06}
Due to the anisotropic interactions between dipoles, dipolar fluids
can present anomalous phase behavior.  Examples of condensed-phase
dipolar systems include ferrofluids, electro-rheological fluids, and
even biological membranes.  Computer simulations have provided useful
information on the structural features and phase transition of the
dipolar fluids. Simulation results indicate that at low densities,
these fluids spontaneously organize into head-to-tail dipolar
``chains''.\cite{Teixeira00,Holm05} At low temperatures, these chains
and rings prevent the occurrence of a liquid-gas phase transition.
However, Tlusty and Safran showed that there is a defect-induced phase
separation into a low-density ``chain'' phase and a higher density
Y-defect phase.\cite{Tlusty00} Recently, inspired by experimental
studies on monolayers of dipolar fluids, theoretical models using
two-dimensional dipolar soft spheres have appeared in the literature.
Tavares {\it et al.} tested their theory for chain and ring length
distributions in two dimensions and carried out Monte Carlo
simulations in the low-density phase.\cite{Tavares02} Duncan and Camp
performed dynamical simulations on two-dimensional dipolar fluids to
study transport and orientational dynamics in these
systems.\cite{Duncan04} They have recently revisited two-dimensional
systems to study the kinetic conditions for the defect-induced
condensation into the Y-defect phase.\cite{Duncan06}

Although they are not traditionally classified as 2-dimensional
dipolar fluids, hydrated lipids aggregate spontaneously to form
bilayers which exhibit a variety of phases depending on their
temperatures and compositions.  At high temperatures, the fluid
($L_{\alpha}$) phase of Phosphatidylcholine (PC) lipids closely
resembles a dipolar fluid.  However, at lower temperatures, packing of
the molecules becomes important, and the translational freedom of
lipid molecules is thought to be substantially restricted.  A
corrugated or ``rippled'' phase ($P_{\beta'}$) appears as an
intermediate phase between the gel ($L_\beta$) and fluid
($L_{\alpha}$) phases for relatively pure phosphatidylcholine (PC)
bilayers.  The $P_{\beta'}$ phase has attracted substantial
experimental interest over the past 30 years. Most structural
information of the ripple phase has been obtained by the X-ray
diffraction~\cite{Sun96,Katsaras00} and freeze-fracture electron
microscopy (FFEM).~\cite{Copeland80,Meyer96} Recently, Kaasgaard {\it
et al.} used atomic force microscopy (AFM) to observe ripple phase
morphology in bilayers supported on mica.~\cite{Kaasgaard03} The
experimental results provide strong support for a 2-dimensional
hexagonal packing lattice of the lipid molecules within the ripple
phase.  This is a notable change from the observed lipid packing
within the gel phase.~\cite{Cevc87}

Although the results of dipolar fluid simulations can not be directly
mapped onto the phases of lipid bilayers, the rich behaviors exhibited
by simple dipolar models can give us some insight into the corrugation
phenomenon of the $P_{\beta'}$ phase.  There have been a number of
theoretical approaches (and some heroic simulations) undertaken to try
to explain this phase, but to date, none have looked specifically at
the contribution of the dipolar character of the lipid head groups
towards this corrugation.  Before we present our simple model, we will
briefly survey the previous theoretical work on this topic.

The theoretical models that have been put forward to explain the
formation of the $P_{\beta'}$ phase have presented a number of
conflicting but intriguing explanations. Marder {\it et al.} used a
curvature-dependent Landau-de Gennes free-energy functional to predict
a rippled phase.~\cite{Marder84} This model and other related
continuum models predict higher fluidity in convex regions and that
concave portions of the membrane correspond to more solid-like
regions.  Carlson and Sethna used a packing-competition model (in
which head groups and chains have competing packing energetics) to
predict the formation of a ripple-like phase.  Their model predicted
that the high-curvature portions have lower-chain packing and
correspond to more fluid-like regions.  Goldstein and Leibler used a
mean-field approach with a planar model for {\em inter-lamellar}
interactions to predict rippling in multilamellar
phases.~\cite{Goldstein88} McCullough and Scott proposed that the {\em
anisotropy of the nearest-neighbor interactions} coupled to
hydrophobic constraining forces which restrict height differences
between nearest neighbors is the origin of the ripple
phase.~\cite{McCullough90} Lubensky and MacKintosh introduced a Landau
theory for tilt order and curvature of a single membrane and concluded
that {\em coupling of molecular tilt to membrane curvature} is
responsible for the production of ripples.~\cite{Lubensky93} Misbah,
Duplat and Houchmandzadeh proposed that {\em inter-layer dipolar
interactions} can lead to ripple instabilities.~\cite{Misbah98}
Heimburg presented a {\em coexistence model} for ripple formation in
which he postulates that fluid-phase line defects cause sharp
curvature between relatively flat gel-phase regions.~\cite{Heimburg00}
Kubica has suggested that a lattice model of polar head groups could
be valuable in trying to understand bilayer phase
formation.~\cite{Kubica02} Bannerjee used Monte Carlo simulations of
lamellar stacks of hexagonal lattices to show that large headgroups
and molecular tilt with respect to the membrane normal vector can
cause bulk rippling.~\cite{Bannerjee02}

Large-scale molecular dynamics simulations have also been performed on
rippled phases using united atom as well as molecular scale
models. De~Vries {\it et al.} studied the structure of lecithin ripple
phases via molecular dynamics and their simulations seem to support
the coexistence models (i.e. fluid-like chain dynamics was observed in
the kink regions).~\cite{deVries05} A similar coarse-grained approach
has been used to study the line tension of bilayer
edges.\cite{Jiang04,deJoannis06} Ayton and Voth have found significant
undulations in zero-surface-tension states of membranes simulated via
dissipative particle dynamics, but their results are consistent with
purely thermal undulations.~\cite{Ayton02} Brannigan, Tamboli and
Brown have used a molecular scale model to elucidate the role of
molecular shape on membrane phase behavior and
elasticity.~\cite{Brannigan04b} They have also observed a buckled
hexatic phase with strong tail and moderate alignment
attractions.~\cite{Brannigan04a} 

The problem with using atomistic and even coarse-grained approaches to
study this phenomenon is that only a relatively small number of
periods of the corrugation (i.e. one or two) can be realistically
simulated given current technology.  Also, simulations of lipid
bilayers are traditionally carried out with periodic boundary
conditions in two or three dimensions and these have the potential to
enhance the periodicity of the system at that wavelength.  To avoid
this pitfall, we are using a model which allows us to have
sufficiently large systems so that we are not causing artificial
corrugation through the use of periodic boundary conditions.

At the other extreme in density from the traditional simulations of
dipolar fluids is the behavior of dipoles locked on regular lattices.
Ferroelectric states (with long-range dipolar order) can be observed
in dipolar systems with non-hexagonal packings.  However, {\em
hexagonally}-packed 2-D dipolar systems are inherently frustrated and
one would expect a dipolar-disordered phase to be the lowest free
energy configuration.  Therefore, it would seem unlikely that a
frustrated lattice in a dipolar-disordered state could exhibit the
long-range periodicity in the range of 100-600 \AA (as exhibited in
the ripple phases studied by Kaasgard {\it et
al.}).~\cite{Kaasgaard03}

Is there an intermediate model between the low-density dipolar fluids
and the rigid lattice models which has the potential to exhibit the
corrugation phenomenon of the $P_{\beta'}$ phase?  What we present
here is an attempt to find a simple dipolar model which will exhibit
this behavior.  We are using a modified XYZ lattice model; details of
the model can be found in section
\ref{sec:model}, results of Monte Carlo simulations using this model
are presented in section
\ref{sec:results}, and section \ref{sec:discussion} contains our conclusions.

\section{2-D Dipolar Membrane}
\label{sec:model}

The point of developing this model was to arrive at the simplest
possible theoretical model which could exhibit spontaneous corrugation
of a two-dimensional dipolar medium.  Since molecules in the ripple
phase have limited translational freedom, we have chosen a lattice to
support the dipoles in the x-y plane.  The lattice may be either
hexagonal (lattice constants $a/b = \sqrt{3}$) or non-hexagonal.
However, each dipole has 3 degrees of freedom.  They may move freely
{\em out} of the x-y plane (along the $z$ axis), and they have
complete orientational freedom ($0 <= \theta <= \pi$, $0 <= \phi < 2
\pi$).  This is essentially a modified X-Y-Z model with translational
freedom along the z-axis.

The potential energy of the system,
\begin{equation}
V = \sum_i \left( \sum_{j \in NN_i}^6
\frac{k_r}{2}\left( r_{ij}-\sigma \right)^2  +  \sum_{j>i} \frac{|\mu|^2}{4\pi \epsilon_0 r_{ij}^3} \left[
{\mathbf{\hat u}_i} \cdot {\mathbf{\hat u}_j} -
3({\mathbf{\hat u}_i} \cdot {\mathbf{\hat
r}_{ij}})({\mathbf{\hat u}_j} \cdot {\mathbf{\hat r}_{ij}})\right]
\right)
\label{eq:pot}
\end{equation}

In this potential, $\mathbf{\hat u}_i$ is the unit vector pointing
along dipole $i$ and $\mathbf{\hat r}_{ij}$ is the unit vector
pointing along the inter-dipole vector $\mathbf{r}_{ij}$.  The entire
potential is governed by three parameters, the dipolar strength
($\mu$), the harmonic spring constant ($k_r$) and the preferred
intermolecular spacing ($\sigma$).  In practice, we set the value of
$\sigma$ to the average inter-molecular spacing from the planar
lattice, yielding a potential model that has only two parameters for a
particular choice of lattice constants $a$ (along the $x$-axis) and
$b$ (along the $y$-axis).  We also define a set of reduced parameters
based on the length scale ($\sigma$) and the energy of the harmonic
potential at a deformation of 2 $\sigma$ ($\epsilon = k_r \sigma^2 /
2$).  Using these two constants, we perform our calculations using
reduced distances, ($r^{*} = r / \sigma$), temperatures ($T^{*} = 2
k_B T / k_r \sigma^2$), densities ($\rho^{*} = N \sigma^2 / L_x L_y$),
and dipole moments ($\mu^{*} = \mu / \sqrt{4 \pi \epsilon_0 \sigma^5
k_r / 2}$). 

To investigate the phase behavior of this model, we have performed a
series of Metropolis Monte Carlo simulations of moderately-sized (34.3
$\sigma$ on a side) patches of membrane hosted on both hexagonal
($\gamma = a/b = \sqrt{3}$) and non-hexagonal ($\gamma \neq \sqrt{3}$)
lattices.  The linear extent of one edge of the monolayer was $20 a$
and the system was kept roughly square. The average distance that
coplanar dipoles were positioned from their six nearest neighbors was
1 $\sigma$ (on both hexagonal and non-hexagonal lattices).  Typical
system sizes were 1360 dipoles for the hexagonal lattices and 840-2800
dipoles for the non-hexagonal lattices.  Periodic boundary conditions
were used, and the cutoff for the dipole-dipole interaction was set to
4.3 $\sigma$.  All parameters ($T^{*}$, $\mu^{*}$, and $\gamma$) were
varied systematically to study the effects of these parameters on the
formation of ripple-like phases.

\section{Results and Analysis}
\label{sec:results}
 
\subsection{Dipolar Ordering and Coexistence Temperatures}
The principal method for observing the orientational ordering
transition in dipolar systems is the $P_2$ order parameter (defined as
$1.5 \times \lambda_{max}$, where $\lambda_{max}$ is the largest
eigenvalue of the matrix,
\begin{equation}
{\mathsf{S}} = \frac{1}{N} \sum_i \left(
\begin{array}{ccc}
        u^{x}_i u^{x}_i-\frac{1}{3} & u^{x}_i u^{y}_i & u^{x}_i u^{z}_i \\
        u^{y}_i u^{x}_i  & u^{y}_i u^{y}_i -\frac{1}{3} & u^{y}_i u^{z}_i \\
        u^{z}_i u^{x}_i & u^{z}_i u^{y}_i  & u^{z}_i u^{z}_i -\frac{1}{3} 
\end{array} \right).
\label{eq:opmatrix}
\end{equation}
Here $u^{\alpha}_i$ is the $\alpha=x,y,z$ component of the unit vector
for dipole $i$.  $P_2$ will be $1.0$ for a perfectly-ordered system
and near $0$ for a randomized system.  Note that this order parameter
is {\em not} equal to the polarization of the system.  For example,
the polarization of the perfect antiferroelectric system is $0$, but
$P_2$ for an antiferroelectric system is $1$.  The eigenvector of
$\mathsf{S}$ corresponding to the largest eigenvalue is familiar as
the director axis, which can be used to determine a privileged dipolar
axis for dipole-ordered systems.  The top panel in Fig. \ref{phase}
shows the values of $P_2$ as a function of temperature for both
hexagonal ($\gamma = 1.732$) and non-hexagonal ($\gamma=1.875$)
lattices.

\begin{figure}[ht]
\centering
\caption{Top panel: The $P_2$ dipolar order parameter as a function of
temperature for both hexagonal ($\gamma = 1.732$) and non-hexagonal
($\gamma = 1.875$) lattices.  Bottom Panel: The phase diagram for the
dipolar membrane model.  The line denotes the division between the
dipolar ordered (antiferroelectric) and disordered phases.  An
enlarged view near the hexagonal lattice is shown inset.}
\includegraphics[width=\linewidth]{phase.pdf}
\label{phase}
\end{figure}

There is a clear order-disorder transition in evidence from this data.
Both the hexagonal and non-hexagonal lattices have dipolar-ordered
low-temperature phases, and orientationally-disordered high
temperature phases.  The coexistence temperature for the hexagonal
lattice is significantly lower than for the non-hexagonal lattices,
and the bulk polarization is approximately $0$ for both dipolar
ordered and disordered phases.  This gives strong evidence that the
dipolar ordered phase is antiferroelectric.  We have repeated the
Monte Carlo simulations over a wide range of lattice ratios ($\gamma$)
to generate a dipolar order/disorder phase diagram.  The bottom panel
in Fig. \ref{phase} shows that the hexagonal lattice is a
low-temperature cusp in the $T^{*}-\gamma$ phase diagram.

This phase diagram is remarkable in that it shows an antiferroelectric
phase near $\gamma=1.732$ where one would expect lattice frustration
to result in disordered phases at all temperatures.  Observations of
the configurations in this phase show clearly that the system has
accomplished dipolar orderering by forming large ripple-like
structures.  We have observed antiferroelectric ordering in all three
of the equivalent directions on the hexagonal lattice, and the dipoles
have been observed to organize perpendicular to the membrane normal
(in the plane of the membrane).  It is particularly interesting to
note that the ripple-like structures have also been observed to
propagate in the three equivalent directions on the lattice, but the
{\em direction of ripple propagation is always perpendicular to the
dipole director axis}.  A snapshot of a typical antiferroelectric
rippled structure is shown in Fig. \ref{fig:snapshot}.

\begin{figure}[ht]
\centering
\caption{Top and Side views of a representative configuration for the
dipolar ordered phase supported on the hexagonal lattice. Note the
antiferroelectric ordering and the long wavelength buckling of the
membrane.  Dipolar ordering has been observed in all three equivalent
directions on the hexagonal lattice, and the ripple direction is
always perpendicular to the director axis for the dipoles.}
\includegraphics[width=5.5in]{snapshot.pdf}
\label{fig:snapshot}
\end{figure}

\subsection{Discriminating Ripples from Thermal Undulations}

In order to be sure that the structures we have observed are actually
a rippled phase and not simply thermal undulations, we have computed
the undulation spectrum,
\begin{equation}
h(\vec{q}) = A^{-1/2} \int d\vec{r}
h(\vec{r}) e^{-i \vec{q}\cdot\vec{r}}
\end{equation}
where $h(\vec{r})$ is the height of the membrane at location $\vec{r}
= (x,y)$.~\cite{Safran94} In simple (and more complicated) elastic
continuum models, Brannigan {\it et al.} have shown that in the $NVT$
ensemble, the absolute value of the undulation spectrum can be
written,
\begin{equation}
\langle | h(q)|^2 \rangle_{NVT} = \frac{k_B T}{k_c |\vec{q}|^4 +
\tilde{\gamma}|\vec{q}|^2},
\label{eq:fit}
\end{equation}
where $k_c$ is the bending modulus for the membrane, and
$\tilde{\gamma}$ is the mechanical surface
tension.~\cite{Brannigan04b}

The undulation spectrum is computed by superimposing a rectangular
grid on top of the membrane, and by assigning height ($h(\vec{r})$)
values to the grid from the average of all dipoles that fall within a
given $\vec{r}+d\vec{r}$ grid area.  Empty grid pixels are assigned
height values by interpolation from the nearest neighbor pixels.  A
standard 2-d Fourier transform is then used to obtain $\langle |
h(q)|^2 \rangle$.

The systems studied in this paper have relatively small bending moduli
($k_c$) and relatively large mechanical surface tensions
($\tilde{\gamma}$).  In practice, the best fits to our undulation
spectra are obtained by approximating the value of $k_c$ to 0.  In
Fig. \ref{fig:fit} we show typical undulation spectra for two
different regions of the phase diagram along with their fits from the
Landau free energy approach (Eq. \ref{eq:fit}).  In the
high-temperature disordered phase, the Landau fits can be nearly
perfect, and from these fits we can estimate the bending modulus and
the mechanical surface tension.

For the dipolar-ordered hexagonal lattice near the coexistence
temperature, however, we observe long wavelength undulations that are
far outliers to the fits.  That is, the Landau free energy fits are
well within error bars for all other points, but can be off by {\em
orders of magnitude} for a few low frequency components.

We interpret these outliers as evidence that these low frequency modes
are {\em non-thermal undulations}.  We take this as evidence that we
are actually seeing a rippled phase developing in this model system.

\begin{figure}[ht]
\centering
\caption{Evidence that the observed ripples are {\em not} thermal
undulations is obtained from the 2-d fourier transform $\langle
|h^{*}(\vec{q})|^2 \rangle$ of the height profile ($\langle h^{*}(x,y)
\rangle$). Rippled samples show low-wavelength peaks that are
outliers on the Landau free energy fits.  Samples exhibiting only
thermal undulations fit Eq. \ref{eq:fit} remarkably well.}
\includegraphics[width=5.5in]{fit.pdf}
\label{fig:fit}
\end{figure}

\subsection{Effects of Potential Parameters on Amplitude and Wavelength}

We have used two different methods to estimate the amplitude and
periodicity of the ripples.  The first method requires projection of
the ripples onto a one dimensional rippling axis. Since the rippling
is always perpendicular to the dipole director axis, we can define a
ripple vector as follows.  The largest eigenvector, $s_1$, of the
$\mathsf{S}$ matrix in Eq. \ref{eq:opmatrix} is projected onto a
planar director axis,
\begin{equation}
\vec{d} = \left(\begin{array}{c}
\vec{s}_1 \cdot \hat{i} \\
\vec{s}_1 \cdot \hat{j} \\
0
\end{array} \right).
\end{equation}
($\hat{i}$, $\hat{j}$ and $\hat{k}$ are unit vectors along the $x$,
$y$, and $z$ axes, respectively.)  The rippling axis is in the plane of
the membrane and is perpendicular to the planar director axis,
\begin{equation}
\vec{q}_{\mathrm{rip}} = \vec{d} \times \hat{k}
\end{equation}
We can then find the height profile of the membrane along the ripple
axis by projecting heights of the dipoles to obtain a one-dimensional
height profile, $h(q_{\mathrm{rip}})$. Ripple wavelengths can be
estimated from the largest non-thermal low-frequency component in the
fourier transform of $h(q_{\mathrm{rip}})$.  Amplitudes can be
estimated by measuring peak-to-trough distances in
$h(q_{\mathrm{rip}})$ itself.

\begin{figure}[ht]
\centering
\caption{Contours of the height-dipole correlation function as a function 
of the dot product between the dipole ($\hat{\mu}$) and inter-dipole 
separation vector ($\hat{r}$) and the distance ($r$) between the dipoles.
Perfect height correlation (contours approaching 1) are present in the 
ordered phase when the two dipoles are in the same head-to-tail line.
Anti-correlation (contours below 0) is only seen when the inter-dipole
vector is perpendicular to the dipoles. }
\includegraphics[width=\linewidth]{height-dipole-correlation.pdf}
\label{fig:CrossCorrelation}
\end{figure}

A second, more accurate, and simpler method for estimating ripple
shape is to extract the wavelength and height information directly
from the largest non-thermal peak in the undulation spectrum.  For
large-amplitude ripples, the two methods give similar results.  The
one-dimensional projection method is more prone to noise (particularly
in the amplitude estimates for the non-hexagonal lattices).  We report
amplitudes and wavelengths taken directly from the undulation spectrum
below.

In the hexagonal lattice ($\gamma = \sqrt{3}$), the rippling is
observed for temperatures ($T^{*}$) from $61-122$.  The wavelength of
the ripples is remarkably stable at 21.4~$\sigma$ for all but the
temperatures closest to the order-disorder transition.  At $T^{*} =
122$, the wavelength drops to 17.1~$\sigma$.

The dependence of the amplitude on temperature is shown in the top
panel of Fig. \ref{fig:Amplitude}.  The rippled structures shrink
smoothly as the temperature rises towards the order-disorder
transition.  The wavelengths and amplitudes we observe are
surprisingly close to the $\Lambda / 2$ phase observed by Kaasgaard
{\it et al.} in their work on PC-based lipids.\cite{Kaasgaard03}
However, this is coincidental agreement based on a choice of 7~\AA~as
the mean spacing between lipids. 

\begin{figure}[ht]
\centering
\caption{a) The amplitude $A^{*}$ of the ripples vs. temperature for a 
hexagonal lattice. b) The amplitude $A^{*}$ of the ripples vs. dipole
strength ($\mu^{*}$) for both the hexagonal lattice (circles) and
non-hexagonal lattice (squares).  The reduced temperatures were kept
fixed at $T^{*} = 94$ for the hexagonal lattice and $T^{*} = 106$ for
the non-hexagonal lattice (approximately 2/3 of the order-disorder
transition temperature for each lattice).}
\includegraphics[width=\linewidth]{properties_sq.pdf}
\label{fig:Amplitude}
\end{figure}

The ripples can be made to disappear by increasing the internal
surface tension (i.e. by increasing $k_r$ or equivalently, reducing
the dipole moment).  The amplitude of the ripples depends critically
on the strength of the dipole moments ($\mu^{*}$) in Eq. \ref{eq:pot}.
If the dipoles are weakened substantially (below $\mu^{*}$ = 20) at a
fixed temperature of 94, the membrane loses dipolar ordering
and the ripple structures. The ripples reach a peak amplitude of
3.7~$\sigma$ at a dipolar strength of 25.  We show the dependence
of ripple amplitude on the dipolar strength in
Fig. \ref{fig:Amplitude}.

\subsection{Non-hexagonal lattices}

We have also investigated the effect of the lattice geometry by
changing the ratio of lattice constants ($\gamma$) while keeping the
average nearest-neighbor spacing constant. The antiferroelectric state
is accessible for all $\gamma$ values we have used, although the
distorted hexagonal lattices prefer a particular director axis due to
the anisotropy of the lattice.

Our observation of rippling behavior was not limited to the hexagonal
lattices.  In non-hexagonal lattices the antiferroelectric phase can
develop nearly instantaneously in the Monte Carlo simulations, and
these dipolar-ordered phases tend to be remarkably flat.  Whenever
rippling has been observed in these non-hexagonal lattices
(e.g. $\gamma = 1.875$), we see relatively short ripple wavelengths
(14 $\sigma$) and amplitudes of 2.4~$\sigma$.  These ripples are
weakly dependent on dipolar strength (see Fig. \ref{fig:Amplitude}),
although below a dipolar strength of $\mu^{*} = 20$, the membrane
loses dipolar ordering and displays only thermal undulations.

The ripple phase does {\em not} appear at all values of $\gamma$.  We
have only observed non-thermal undulations in the range $1.625 <
\gamma < 1.875$.  Outside this range, the order-disorder transition in
the dipoles remains, but the ordered dipolar phase has only thermal
undulations.  This is one of our strongest pieces of evidence that
rippling is a symmetry-breaking phenomenon for hexagonal and
nearly-hexagonal lattices.

\subsection{Effects of System Size}
To evaluate the effect of finite system size, we have performed a
series of simulations on the hexagonal lattice at a reduced
temperature of 122, which is just below the order-disorder transition
temperature ($T^{*} = 139$).  These conditions are in the
dipole-ordered and rippled portion of the phase diagram.  These are
also the conditions that should be most susceptible to system size
effects.

\begin{figure}[ht]
\centering
\caption{The ripple wavelength (top) and amplitude (bottom) as a
function of system size for a hexagonal lattice ($\gamma=1.732$) at $T^{*} =
122$.}
\includegraphics[width=\linewidth]{SystemSize.pdf}
\label{fig:systemsize}
\end{figure}

There is substantial dependence on system size for small (less than
29~$\sigma$) periodic boxes.  Notably, there are resonances apparent
in the ripple amplitudes at box lengths of 17.3 and 29.5 $\sigma$.
For larger systems, the behavior of the ripples appears to have
stabilized and is on a trend to slightly smaller amplitudes (and
slightly longer wavelengths) than were observed from the 34.3 $\sigma$
box sizes that were used for most of the calculations.

It is interesting to note that system sizes which are multiples of the
default ripple wavelength can enhance the amplitude of the observed
ripples, but appears to have only a minor effect on the observed
wavelength.  It would, of course, be better to use system sizes that
were many multiples of the ripple wavelength to be sure that the
periodic box is not driving the phenomenon, but at the largest system
size studied (70 $\sigma$ $\times$ 70 $\sigma$), the number of dipoles
(5440) made long Monte Carlo simulations prohibitively expensive. 

\section{Discussion}
\label{sec:discussion} 

We have been able to show that a simple dipolar lattice model which
contains only molecular packing (from the lattice), anisotropy (in the
form of electrostatic dipoles) and a weak surface tension (in the form
of a nearest-neighbor harmonic potential) is capable of exhibiting
stable long-wavelength non-thermal surface corrugations.  The best
explanation for this behavior is that the ability of the dipoles to
translate out of the plane of the membrane is enough to break the
symmetry of the hexagonal lattice and allow the energetic benefit from
the formation of a bulk antiferroelectric phase.  Were the weak
surface tension absent from our model, it would be possible for the
entire lattice to ``tilt'' using $z$-translation.  Tilting the lattice
in this way would yield an effectively non-hexagonal lattice which
would avoid dipolar frustration altogether.  With the surface tension
in place, bulk tilt causes a large strain, and the simplest way to
release this strain is along line defects.  Line defects will result
in rippled or sawtooth patterns in the membrane, and allow small
``stripes'' of membrane to form antiferroelectric regions that are
tilted relative to the averaged membrane normal.

Although the dipole-dipole interaction is the major driving force for
the long range orientational ordered state, the formation of the
stable, smooth ripples is a result of the competition between the
surface tension and the dipole-dipole interactions.  This statement is
supported by the variation in $\mu^{*}$.  Substantially weaker dipoles
relative to the surface tension can cause the corrugated phase to
disappear.

The packing of the dipoles into a nearly-hexagonal lattice is clearly
an important piece of the puzzle.  The dipolar head groups of lipid
molecules are sterically (as well as electrostatically) anisotropic,
and would not be able to pack hexagonally without the steric
interference of adjacent molecular bodies.  Since we only see rippled
phases in the neighborhood of $\gamma=\sqrt{3}$, this implies that
there is a role played by the lipid chains in the organization of the
hexagonally ordered phases which support ripples in realistic lipid
bilayers.

The most important prediction we can make using the results from this
simple model is that if dipolar ordering is driving the surface
corrugation, the wave vectors for the ripples should always found to
be {\it perpendicular} to the dipole director axis.  This prediction
should suggest experimental designs which test whether this is really
true in the phosphatidylcholine $P_{\beta'}$ phases.  The dipole
director axis should also be easily computable for the all-atom and
coarse-grained simulations that have been published in the literature.

Our other observation about the ripple and dipolar directionality is
that the dipole director axis can be found to be parallel to any of
the three equivalent lattice vectors in the hexagonal lattice.
Defects in the ordering of the dipoles can cause the dipole director
(and consequently the surface corrugation) of small regions to be
rotated relative to each other by 120$^{\circ}$.  This is a similar
behavior to the domain rotation seen in the AFM studies of Kaasgaard
{\it et al.}\cite{Kaasgaard03}   

Although our model is simple, it exhibits some rich and unexpected
behaviors.  It would clearly be a closer approximation to the reality
if we allowed greater translational freedom to the dipoles and
replaced the somewhat artificial lattice packing and the harmonic
``surface tension'' with more realistic molecular modeling
potentials.  What we have done is to present an extremely simple model
which exhibits bulk non-thermal corrugation, and our explanation of
this rippling phenomenon will help us design more accurate molecular
models for corrugated membranes and experiments to test whether
rippling is dipole-driven or not.
\clearpage
\bibliography{ripple}
\printfigures
\end{document}
Revision:	3091
Committed:	Tue Dec 26 21:50:50 2006 UTC (18 years, 10 months ago) by xsun
Content type:	application/x-tex
File size:	32281 byte(s)
Log Message:	Making changes to satisfy the reviewers