\documentclass[10pt]{article}


\usepackage{helvet}
\renewcommand{\rmdefault}{phv}

\usepackage{amssymb}        
\usepackage[pctex32]{graphics}


\usepackage{amsmath}
\usepackage{amsfonts}


\setlength{\textwidth}{6.5in}
\setlength{\textheight}{9.0in}

%\setlength{\oddsidemargin}{-0.25in}
\setlength{\oddsidemargin}{0in}
%\setlength{\evensidemargin}{0.25in}
\setlength{\topmargin}{-.50in}





\renewcommand{\topfraction}{0.85}
\renewcommand{\textfraction}{0.1}
\renewcommand{\floatpagefraction}{0.75}








\renewcommand{\baselinestretch}{1.0}
\def\ds{\displaystyle}
\newcommand\bsig{\mbox{\boldmath${\sigma}$}}
\newcommand\beps{\mbox{\boldmath${\epsilon}$}}
\newcommand\bnab{\mbox{\boldmath${\nabla}$}}
\newcommand\bcdot{\mbox{\boldmath${\cdot}$}}
\newcommand\btimes{\mbox{\boldmath${\times}$}}
\newcommand\bxi{\mbox{\boldmath${\xi}$}}


\newtheorem{Thm}{Theorem}[section]


\newcommand{\supp}[0]{\mathrm{supp}}



\def\e{\epsilon}
\def\s{\sigma^*}
\def\D{\partial}
\def\->{\rightarrow}
\def\dom{{\cal D}_h}
\def\O{\Omega}
\def\o{\omega}
\def\t{\tau}
\def\Z{\cal{Z}}
\def\S{{\cal S}_\e}
\def\F{\cal{F}}



\def\R{{\mathbb R}}




\newcommand{\vv}{\textbf{v}}
\newcommand{\xx}{\textbf{x}}
\newcommand{\ff}{\textbf{f}}
\newcommand{\KK}{\textbf{K}}
\newcommand{\VV}{\textbf{V}}
\newcommand{\kk}{\textbf{k}}
\newcommand{\II}{\textbf{I}}




\begin{document}


% Tony - what will knowing more about dependence of
%EM props on microstructure how will that help worbot,
%remote sensing of thickness what other antarctic science
%studies will we be helping
%
%Maloclm - what potential is there for leaving 
%cross-borehole arrays in the ice -- in Septemeber-October
%early spring time frame and transmitting data to
%satellites?


%% broader impact --- ETOPIM
% phys. math geophysics, bio, EE, etc.

%% broader impact book - incorporate polar science
%% examples into the calculus curriculum!
%% everyone loves penguins - calcu students don's expect
%% to see penguins in calc class, AND be shown how their
%% math is related to studying the biology of the algae
%% living in the sea ice 

%webcast back to calc classes at U of U and UAF 



\pagestyle{empty}


\begin{center}

%{\large \bf PROJECT SUMMARY \\ 
%\vspace{1ex}
%Collaborative Research: Probing Antarctic Sea Ice  \\Processes with %Low-Frequency Electromagnetic Waves}






%{\large \bf PROJECT SUMMARY \\ 
%\vspace{-3ex}
%Collaborative Research: 
%Probing Internal Structure and Transport \\ in Sea ice 
%to Monitor Climate Change}


%{\large \bf PROJECT SUMMARY \\ 
%\vspace{1ex}
%Collaborative Research: 
%Climate Change and \\ Electromagnetic Monitoring of
%Transport in Sea Ice}




{\bf {\Large PROJECT SUMMARY}}\\


\vspace{2ex}



{\bf {\large Phase Transitions in Composite Media }}\\ 

%\vspace{0.5ex}




\end{center}

%\vspace{-1.2ex}


Transport in disordered composites is a problem which arises throughout the sciences and engineering and has attracted significant theoretical, computational, and experimental interest. One of the key features of these types of problems is the critical dependence of the effective transport properties on system parameters, such as volume fraction, component contrast ratio, applied field strength, etc. In recent years a broad range of mathematical techniques have been developed to study phase transitions exhibited by such composites, revealing features which are virtually ubiquitous in disordered systems.

%Transport in disordered two--component composite media is a problem which arises throughout the sciences and engineering and has attracted significant theoretical, computational, and experimental interest. One of the interesting and key features of these types of problems is the critical dependence of the effective transport properties on system parameters, such as volume fraction, component contrast ratio, applied field strength, temperature, etc. In fact, the effective properties of numerous disordered composite systems exhibit critical behavior as parameters vary. Examples of such disordered composites include porous media, electrorheological (ER) fluids, sea ice, percolating conductors, and many biocomposites. In recent years a broad range of mathematical techniques have been developed to study several classes of phase transitions exhibited by disordered composites. These techniques have revealed many universal features which are virtually ubiquitous in disordered systems, arising in important problems in both pure and applied mathematics, such as disordered quantum/mesoscopic conductors, quantum chaos, neural networks, random graph theory, atomic energy resonances, and even the Riemann hypothesis.

Here we propose to conduct fundamental, mathematical studies of phase transitions and ground state geometry in two phase composite media in conjunction with state of the art numerical simulations and small scale laboratory experiments. In particular, we will utilize and extend techniques of statistical mechanics and random matrix theory to study conduction/insulation phase transitions in classical two component conductors, liquid/solid phase transitions in electrorheological (ER) suspensions, and microstructural and transport transitions in sea ice.
Toward this goal, we present here a statistical mechanics framework 
for transport in binary composites, placing this class of problems 
on an equal mathematical footing with the Ising model. This approach is
developed from first principles in physics and is based on
an analytic representation of the system energy, which follows from a resolvent formula for the electric field involving a sparse, self adjoint random matrix. 
%which separates the constituent parameters from 
%the applied field strength and geometric information.
Key features of the proposed work include:

\begin{itemize}
%\vspace{-1ex}
\item Characterization of the critical electric field 
in the ER liquid/solid phase transition. 
Preliminary analytical, numerical, and experimental results 
indicate that the energy landscape exhibits a hilly terrain 
associated with metastable states.  Critical behavior 
appears to occur when the applied field
is strong enough to overcome
these local energy minima, a phenomenon which apparently goes beyond 
classical phase transitions in statistical mechanics.

\vspace{-1ex}
%
\item Investigation of diffusion limited
aggregation (DLA) clusters as configurations that minimize a
free energy.  We will study the relationship between 
metal sphere ER fractal ground states and DLA clusters 
through comparisons of spectral properties and statistical mechanical 
studies.




\vspace{-1ex}
\item Spectral analysis of connectedness and transport
transitions in random resistor networks, combining techniques 
of random matrix theory and the theory of polynomials orthogonal to the spectral measure, thereby coupling geometric fluctuations with universal eigenvalue fluctuations.

\vspace{-1ex}

\item Investigation of the connectedness transition in the
brine microstructure of sea ice, and related critical behavior 
of transport properties. X-ray computed tomography can provide us
with unprecedented data on its thermal evolution,
which we will examine through the lens of random matrix theory 
and the corresponding evolution of spectral statistics. 
This approach will yield new insights into
fluid and electrical transport in sea ice, which 
mediate a broad range of climatological and biogeochemical
processes, and are central to monitoring the impact of
climate change on the polar marine environment.

\end{itemize}


The {\bf intellectual merit} of our proposal is that it will
allow us to identify and study, with powerful analytic and computational
methods, key features of phase transitions in composite media 
in diverse, applied settings. Our proposed work will bring to bear 
novel approaches to studying such problems, taken from random matrix
theory and statistical mechanics. Through unique outreach
and educational opportunities, as well as the close connection
of our work to topics of broad public interest such as climate
change, we are in a strong position to facilitate the 
{\bf broader impact} of our proposed research. The postdocs and students
will be involved in helping to understand the rich behavior of 
certain types of composites, which are important to a number of
emerging applications in science and technology. Moreover, 
they will be immersed in a highly interdisciplinary research 
environment, spanning analysis, computation, laboratory
experimentation, and even polar field work. This type of project
will attract young researchers across disciplines to modern applied mathematical research, and highlight the importance of mathematics
in such areas as climate research and the development of high-tech
composite systems. 






\eject







\setcounter{page}{1} 
\pagestyle{plain} 



\noindent
{\Large {\bf PROJECT DESCRIPTION}}





\vspace{2ex}

\noindent
{\large {\bf Results of Prior NSF Support}}

\vspace{1ex}




\noindent {\bf AWARD:} NSF Division of Atmospheric Sciences
Grant DAS-0222171, 7/02 - 9/06; 
Collaboration in Mathematical Geosciences (CMG): 
{\it Microstructural Controls on Transport Processes 
in Geophysical Systems.} PI: $\;$K. M. Golden, U. of Utah (UU),  
co-PI: $\;$H. Eicken, U. of Alaska, Fairbanks (UAF), \$733,000, 
including a Research Experiences for Undergraduates (REU) 
Supplement to UU. %(\$63,000).

\vspace{0.5ex}

\noindent 
{\bf SUMMARY OF RESULTS:} 
The vertical fluid permeability $k$ of sea ice constrains a broad 
range of processes, such as the growth and decay of seasonal ice, 
the evolution of summer ice albedo, and biomass build-up. 
%Such processes are critical to how sea ice and associated ecosystems 
%respond to climate change. 
However, studies of $k$ and its dependence 
on brine porosity $\phi$ and microstructure were sparse. 
We developed a unified theory for $k(\phi)$ which closely captures 
lab and field data \cite{Golden:2007:GRL,Golden:2006:MM,Zhu:2006:AG}. 
Our analytical and numerical modeling is based on rigorous bounds 
\cite{Golden:2006:MM}, lattice and continuum percolation theory 
\cite{Golden:S-2238,Golden:PRL-2923,Golden:2007:GRL}, 
hierarchical models \cite{Golden:2007:GRL}, and
random pipe networks \cite{Zhu:2006:AG}. 
We found that sea ice displays universal transport properties 
remarkably similar to crustal rocks. 
X-ray computed tomography 
and pore structure analysis provided an unprecedented look at the 
brine phase and its connectivity \cite{Golden:2007:GRL}.
These advances have also
allowed us to predict anisotropic percolation thresholds in sea ice,
and to compute, for the first time in natural materials, various
connectedness functions of percolation theory
\cite{Pringle:JGR:2009}. 
In related work, we obtained bounds 
on the complex permittivity of sea ice at 50 MHz, and inverted the 
data to get bounds on brine porosity \cite{Gully:2007:Physica}.
%a first step toward the proposed work here. 
Finally, we used percolation models to help
explain fluid production and anomalies in resistivity data 
in rocks deep in the mid-crust 
of the Southern Alps of New Zealand \cite{Vry:NG:2009}.

%Our results permit more realistic 
%representations of sea ice in global climate and biogeochemical models. 


Our work on permeability has had broad impact.
From 2007 through November 2009, Golden has been invited to
give over 45 lectures on this work at universities
and major conferences in the
US and abroad. This work has also received extensive media
coverage, including {\it Geotimes}, {\it Ars Technica},  
{\it Electrical Engineering Times}, a two page profile
of Golden 
in {\it Science} (3 April 2009), 
radio, television and web interviews, and numerous invited lectures 
to general audiences and student groups. For example, 
Golden was invited by 
the AMS to give a Congressional Luncheon Briefing on Capitol Hill 
in Washington D. C. in November 2007.
%(www.ams.org/government/congress-briefing-nov07.html).


\vspace{0.5ex}

\noindent 
{\bf CONTRIBUTION TO HUMAN RESOURCES:} This award has provided 
partial support for three graduate students, J. Miner and 
L. Backstrom at UAF, and B. Murphy at UU, a post-doc, 
D. Pringle at UAF, two programmers, and 
six UU undergraduates from a range of majors. Support 
included travel between UU and UAF, 
and Arctic field experience for five of the undergraduates. 
All students were involved in publications and giving 
lectures. For example, Amy Heaton
\cite{Golden:2007:GRL,Golden:2006:MM} began with
Golden in his freshman Calculus class (as did many of the other
REU students), and
has given presentations 
%at an Undergraduate Math Res. Conf. at 
at Ohio State (2003), 
in Physics at Victoria U., NZ (2004), 
the Utah State Legislature (2004), and with Golden 
in the US Congress in June 2003. She completed
a Ph.D. in theoretical quantum chemistry from UU in 2009. 
%(www.ams.org/government/cnsfex03.html). 
Ali Jabini went on to complete an M.S. in EE 
at Columbia, and Troy Finlayson completed an M.D. at UU.
In 2007 
Kellen Petersen won the Research Scholar Award, the highest 
undergraduate research honor in the 
College of Science at UU, and is now in the Math Ph.D.
program at Courant Institute. In 2008 Adam Gully won the same award. 
%Four of these students, and many other REU's, took Calculus with Golden.











\vspace{1ex}



\noindent {\bf AWARD:} NSF Division of Mathematical Sciences
Grant DMS-0537015, 9/1/05 - 8/31/09 (and a one year
NSF granted extension), UU,
{\it Analysis and Computation of Electromagnetic 
Transport in Composite Materials}.
PI: $\;$K. M. Golden,  co-PI: $\;$D. Dobson,
\$403,000; including an REU Supplement. % (\$61,000).

\vspace{0.5ex}

\noindent 
{\bf SUMMARY OF RESULTS:} We have been investigating novel 
approaches to some challenging issues in electromagnetic 
transport in composites. %Several results have 
%been obtained in the effort to extend the analytic 
%continuation method for bounding the effective complex permittivity 
%beyond quasistatics to dynamics \cite{Simeonova:2009,Simeonova:2009B}.  
We've obtained rigorous and numerical results on a new
effective complex permittivity $\epsilon^*$ 
defined for finite frequency, and how it behaves 
in certain limits \cite{Simeonova:2009,Simeonova:2009B}. 
Former graduate student Lyubima Simeonova studied effective medium theory
for wave propagation in random cell materials, outside of the usual
long wavelength limit.  Working with a natural
extension of the quasistatic effective medium parameter, which is no
longer necessarily constant, she used probability arguments along
with PDE analysis to provide estimates of $\epsilon^*$.
She conducted
numerical experiments illustrating the limits of the bounds,
and developed an optimization scheme 
yielding the probability density
for the class of materials which maximizes the spatial 
average of $\epsilon^*$ \cite{Simeonova:2009C,Simeonova:2009B}.  
The investigation
into optimizing effective material properties led to 
the construction of a periodic composite with
properties similar to that of a slab of negative refractive index
material \cite{Dobson:2009:AMO}. 
Simeonova now holds a postdoctoral position at Texas A\&M.




Electrorheological (ER) fluids are particle suspensions that
exhibit a rapid fluid/solid transition when an applied electric
field exceeds a critical value. Particles can form columns
%with periodic microstructure
or fractal networks.
Graduate student Hwanyong Lee developed a computational
method for calculating the free energy for 
``two-dimensional'' ER fluids with disc-shaped particles.  
The effective dielectric constant $\epsilon^*$ in the quasistatic
limit can be computed by solving the
electrostatic problem (numerically, following Greengard's
fast multipole method) over the domain and finding the
spatially averaged displacement field, and the free energy can be
computed from $\epsilon^*$. Error estimates were also developed.
%Our approach to solving the
%electrostatic problem follows classical methods of Greengard, using
%integral equations on a large number of disc-shaped particles.
%Error estimates were developed 
%which indicate how accurately
%each electrostatics problem must be solved, and for how many values
%of s one much compute F(s), 
%in order to find the free energy up to a given tolerance. 
This method is awaiting numerical implementation
before publication. Lee is expected to graduate in May 2010.
REU Peter Sommerkorn, who graduated in May 2009 
with majors in Mechanical Engineering and Math,
conducted experiments on ER fluids, imaging 
their microstructural evolution with
cameras and video, and measuring shear viscosity in
parallel plate and circular geometries. Peter led
5 other ME seniors in these experiments for their
Senior Design projects, and they won the
2009 Schlumberger Mechanical Engineering Project Award.


Graduate student Russell Richins investigated
applications of a new minimization-based
formulation of time-harmonic wave equations (due to Graeme Milton and
co-authors).  He has developed and implemented a finite element
method based on this strategy, complete with error estimates.
It is anticipated that this minimization principle will yield another
approach to obtaining bounds on effective material parameters for wave
propagation, at least in the case of dissipative materials. Richins
is expected to graduate in May 2010.




Graduate student Ben Murphy's work has been on connecting 
the critical properties of transport in composites to phase transitions 
in statistical mechanics and random matrix theory
\cite{Murphy_Stat_Mech,Murphy_RMT}. This work helps forms the basis
for this proposal and will be described below.
He is expected to graduate in May 2011.


A composite where many key issues of EM 
transport arise is sea ice $-$
pure ice with brine inclusions.
%which shares connectivity 
%threshold behavior with 
%ER fluids. 
In Sept.-Oct. 2007, Golden and math student 
Adam Gully participated in the Sea Ice Physics and 
Ecosystem eXperiment (SIPEX) in Antarctica aboard the 
Australian icebreaker {\it Aurora Australis}. We used an 
earth resistance tester in a Wenner array to obtain DC 
resistivity profiles of the ice, and adapted the set-up 
to measure the vertical conductivity of cores, in the
first direct measurements of this type \cite{Golden:2010:DSR}. 
We found strong evidence of a percolation threshold 
for electrical conductivity \cite{Golden:GRL:2010},
with implications for monitoring sea ice thickness.
We developed percolation, hierarchical and network
models \cite{Zhu:PB:2010} to explain our results
\cite{Golden:2010:DSR,Golden:GRL:2010}, obtaining excellent
agreement with our data. 
We also 
%a wave generator and oscilloscope to 
measured the complex permittivity of a sea ice core from 
10 KHz to 15 MHz. 
We made the first measurements of fluid permeability in 
Antarctic pack ice, 
%which are consistent with our 
%lognormal pipe bounds \cite{Golden:2007:GRL,Golden:2006:MM}, 
and conducted tracer experiments on ice blocks. 
We found a new percolation threshold for fluid flow in 
granular sea ice, 
which has broad implications for climatological and biological
processes \cite{Golden:2010:FTASI}.


The sea ice work, particularly Golden and Gully's 
findings during the SIPEX expedition, is having broad impact. 
The AMS has featured Golden's Antarctic work in a ``Mathematical Moment'' 
for the general public on its website, including an four-part interview,
and Golden was interviewed by satellite phone from Antarctica
for Public Radio during the expedition. 
Golden was invited to write an extensive 
article for the {\it Notices of the AMS}, which was featured
on the cover of the May 2009 issue \cite{Golden:NAMS:2009}. 
Golden has produced a DVD giving an overview of the science 
on SIPEX, and the excitement and beauty of venturing to 
Antarctica. This footage has been shown on TV
and to many student, scientific, and general audiences,
including the US Congress and the largest visual arts
college in the US, Columbia College Chicago.
Golden and Gully were recently featured on a 30 minute
show on KUED Public TV, {\it Utah Conversations with Ted Capener}:
``Sea ice, mathematics and global climate change: Is there a link?
... A conversation with Ken Golden.''
(see www.math.utah.edu/$\sim$golden/media.html for details).



\vspace{0.5ex}

\noindent 
{\bf CONTRIBUTION TO HUMAN RESOURCES:} 
Partial support for 5 math grad students Simeonova, 
Murphy, Richins, Lee, and Gully, 
and 6 undergrads, Olakunle Eso, Megan Morris, Dave Arcilesi, 
Christian Sampson, Peter Sommerkorn, and Kellen Petersen. 
%Simeonova will finish her Ph.D. in spring 2009, and Murphy is expected
%to finish his Ph.D. in spring 2010. 
%Most of these students have participated 
%in a weekly EM seminar with Golden, Dobson, and Zhu, 
%where the students usually give the lectures. 







\vspace{1ex}




\noindent {\bf AWARD:} NSF Div. Math. Sci.
Grant DMS-0602219, 5/17/2006 -- 7/31/2011, UU,
{\it EMSW21-VIGRE: Vertical Integration in Mathematics at the University of Utah}.
PI: $\;$A. Bertram,  coPI: $\;$D. Dobson, Senior Personnel:
     F. Adler, E. Cherkaev, K. M. Golden, N. Korevaar, G. Savin, 
     K. Schmitt, P. Trapa, and J. Zhu, \$3,500,000.
The UU Math. Dept. received two consecutive
VIGRE awards in 2001 and 2006, supporting 
pre-college, undergraduate, graduate, and postdoctoral 
programs. As Director of Undergrad. Studies (2002-8), 
and as REU Coordinator (2003-7), Golden was
responsible for the undergraduate components. 
%Since 2001, they have so far supported 156 undergrads in REU programs,
%58 grad students, and 17 postdocs.
Since 2001, they have so far supported 179 undergrads in REU
programs, 55 grad students, and 19 postdocs.
%Since 2001, they have so far supported 114 undergrads in 
%REU programs, 45 grad students, and 15 postdocs. 
Golden led efforts to broaden our degrees, introducing a new major in 
Applied Math designed to promote double 
majoring and interdisciplinary studies. 
Rapid growth in the number of math majors has resulted.
%Since 2002 the number of math majors (and pre-majors) at UU has roughly doubled. 













\section{Introduction}




\begin{figure}                \centering
\includegraphics{Figure_1.eps}
   \caption{(a) Fractal network of metal spheres in an 
             ER fluid (W. Wen).
            (b) Horizontal cross-section of a brine channel 
             in Arctic sea ice (K. Golden).
             (c) Diffusion limited 
     aggregation cluster \cite{Meakin:1998,Bunde:FDS-91}.}
      
    \label{erfluid_brine_dla}
        
\end{figure}


Composite media arise naturally throughout the physical and 
biological sciences, and are employed in a broad range of
engineering and technological applications. 
The behavior of those media exhibiting a critical transition
as system parameters are varied 
is particularly challenging to describe physically,
and to predict mathematically. 
For example electrorheological (ER) fluids, which are suspensions 
of spherical particles in a dielectric liquid host, 
exhibit a liquid to solid phase transition 
as the strength $E_0$ of an applied electric field  
surpasses a critical value $E_c$
\cite{Halsey:S-761,Tao-EF,Wen:APL-2147,Wen:PF-1826}. 
Plastic or glass dielectric spheres, with diameters
on the scale of tens of microns,
form clusters and then chains which coalesce into columns
with periodic lattice arrangements of the spheres
as the field increases.
Metal spheres form fractal net structures as shown in 
Figure \ref{erfluid_brine_dla} (a). This system undergoes
an electrically-induced transition in the connectedness of the spheres.
With an increase in viscosity 
of these suspensions by several 
orders of magnitude within a few 
milliseconds \cite{Wen_etal}, they have been 
used in clutches, brakes, micro-fluidic valves, human prosthetics, and
active mechanical elements capable of responding to 
environmental variations.

Another composite which displays complex critical behavior is sea ice,
consisting of pure ice with sub-millimeter brine inclusions, 
whose volume fraction $\phi$,
geometry, and connectedness vary significantly with temperature $T$.
The polar ice packs are both indicators and agents of climate change. 
They also host extensive algal and bacterial communities, which live 
in the brine inclusions and sustain life 
in the polar oceans. Fluid flow through porous sea ice mediates 
a broad range of processes such as the growth and decay of 
seasonal ice, the evolution of surface melt ponds and ice pack 
reflectance, and biomass build-up \cite{Golden:NAMS:2009}.
In \cite{Golden:S-2238} we identified the {\it on-off switch}
for fluid transport in sea ice.
For brine volume fractions $\phi$ below about $5\%$, columnar
sea ice is effectively impermeable to fluid flow, while for 
$\phi$ above $5\%$, it is increasingly permeable.
This critical brine volume fraction  $\phi_c \approx 5\%$ 
corresponds to a critical temperature 
$T_c \approx -5^{\circ}$~C for a typical
bulk sea ice salinity of 5 parts per thousand, which is known as the 
{\it rule of fives}. Fluid
flow is facilitated by
{\it brine channels} $-$ connected brine structures
ranging in scale from a few centimeters for horizontal
cross-sections, as shown in Figure \ref{erfluid_brine_dla} (b), 
to a meter or more in the vertical direction. 
The critical behavior of fluid flow through sea ice was 
postulated in \cite{Golden:S-2238} to result from a 
temperature-driven transition in connectedness 
of the brine microstructure, or a {\it percolation threshold}. 
We predicted the 
5\% critical brine volume fraction in sea ice
by noting the close similarity of its microstructure 
to that of stealthy, radar absorbing
materials, and adapted for sea ice a continuum percolation model for
compressed powders \cite{Kusy:N-58,Kusy:JAP-5301}
used in the design of these materials. 

Lattice and continuum percolation models have 
been used to study a 
broad range of materials including rocks 
\cite{Bourbie:JGR-11524,Broadbent:PCPS-629}, 
semiconductors \cite{Efros-84}, thin films 
\cite{Davis:OC-70}, glacial ice \cite{Enting:1985:LSM}, 
polycrystalline metals \cite{Chen:PRL:2007}, 
radar absorbing coatings \cite{Kusy:N-58}, and 
carbon nanotube composites \cite{Kyrylyuk:PNAS:2008}.  
In the simplest case of the two dimensional square lattice
\cite{Stauffer-92,Torquato:RHM-02},
the bonds are 
open with probability $p$ and closed with 
probability $1-p$.
Connected sets of open bonds 
are called open clusters. 
The average cluster size
grows as $p$ increases, and
%The striking feature of the model is that 
there is a critical probability $p_c$, $0<p_c<1$,
called the {\it percolation threshold}, 
where an infinite cluster of open bonds first appears.
In dimension $d=2$, $p_c = \frac{1}{2}$, and in $d=3$,
$p_c \approx 0.25$. Now consider transport through
the associated {\it random resistor network}, where the
bonds are assigned
electrical conductivities $\sigma_1$
with probability $1-p$, and $\sigma_2=1$
with probability $p$. Then the
effective
conductivity $\sigma^*(p)$ exhibits critical behavior
as $h=\sigma_1/\sigma_2 \rightarrow 0$. 
For $p<p_c$,
$\sigma^*(p)=0$, while $\sigma^*(p)>0$ for $p>p_c$, with  
$\sigma^*(p) \sim (p - p_c)^t$, 
as $p \rightarrow p_c^+$, where $t$ is the conductivity
critical exponent, believed to be {\it universal} for lattices,
depending only on dimension.


The critical behavior of the random resistor 
network is reminiscent of a phase transition in statistical
mechanics 
\cite{Clerc:AP-191,Bergman:SSP-147,Efros:PSS-475,Hong:PRB-4564}, 
like that exhibited by an
Ising ferromagnet \cite{Thompson-88,Christensen:2005}
around its Curie point at a critical
temperature $T=T_c$, as the applied magnetic field $H \rightarrow 0$.
The formulations and basic physics of these two 
classes of problems are nevertheless quite different.
In \cite{Golden:JMP-5627,Golden:PRL-3935}, however, it was 
observed that the {\it analytic continuation method} for 
bounding effective transport properties of composites
\cite{Bergman:PRL-1285,Milton:APL-300,Golden:CMP-473}
provides a mathematical link between them, through
the Lee-Yang Theorem in statistical
mechanics \cite{Lee:PR-411,Baker:PRL-990,Baker-1990}.
This theorem, which states that the zeros of the
partition function lie on the unit circle in the activity
variable,
yields a logarithmic integral for the free
energy of an Ising model, and a Stieltjes representation
for the magnetization $M(H)$.
In the analytic continuation method, $m=\sigma^*/\sigma_2$ is 
an analytic function
of $h$, taking the upper half plane to the upper half plane.
A Stieltjes integral representation for 
$F(s)=1-m(h)$, where $s=1/(1-h)$, is used to obtain 
rigorous bounds on $\sigma^*$ given
microstructural information. This formula
is based on a resolvent representation for the 
electric field, involving a self adjoint random matrix
(or operator). All the geometry of the 
composite is incorporated through a spectral measure of this
matrix. Through these {\it mathematical} parallels
the two classes of problems can be placed on an equal footing.

This analytic parallel has been exploited
to show that the critical exponents of transport in both 
lattice and continuum percolation models
obey the classical scaling relations of statistical mechanics
\cite{Golden:PRL-3935}. However, without establishing the physics 
underlying such a connection, or exploring 
further the mathematics common to both types
of problems, the potential of such methods remains largely 
untapped.

{\bf Here we propose} to do just that, namely, to 
build the statistical mechanics framework for studying
phase transitions in composites, and to explore and exploit random 
matrix theory as it pertains to critical behavior.
We will use a few key questions to guide our investigations.
A central, unresolved issue in the theory of ER fluids is how
to characterize the critical electric field required to induce
the transition. Simple considerations suggest that the application
of any electric field could drive the system to a connectedness
transition, but this is not observed in simulations or experiments.
That is, why is there a critical electric field required for the
transition at all? 
The ER fluid is an almost ideal testing ground where a statistical
mechanics approach, adapted to macroscopic composite media, 
can yield fundamental insights into the nature of the transition.
Another related question is the ground state, particularly
for metal spheres. The resulting fractal nets are very similar in structure
to diffusion limited aggregation (DLA) clusters 
\cite{Meakin:1998,Bunde:FDS-91},
and they may share
some deep mathematical connections. A DLA cluster is shown in 
Figure \ref{erfluid_brine_dla} (c). These aggregates
are defined from a purely probabilistic set of rules, yet we will 
explore the possibility of characterizing them as energy minimizing 
ground states, via the connection to the ER fluid.
Our ER fluid studies will be done in close collaboration with
Ping Sheng 
(Dr. William M. W. Mong Chair Professor of Nanoscience)
and Weijia Wen (Associate Professor of Physics)
of the Hong Kong University of Science and Technology.


We also plan to explore the connectedness transition in the brine
microstructure of sea ice, and more generally the percolation 
threshold in random resistor networks through analysis
of associated random matrices and their spectral properties. 
Characterizing temperature-driven microstructural transitions 
in sea ice through the evolution of spectral statistics would
provide fundamental new insights into its transport behavior.
Such information can be obtained from X-ray 
tomography, and would open a new chapter in the application
of random matrix theory to the analysis of geophysical and 
other complex materials. We will investigate in detail
the spectral characteristics of the resistor network, as well as its
description in terms of our statistical mechanics framework, 
as a key model for sea ice microstructure and transport properties.




The recent development of numerical methods 
\cite{Wang2001, Myroshnychenko2005, Cheng2008} to calculate 
the interaction of fields with complex, 3-d geometries of 
particles in a host matrix provides an invaluable tool for 
simulating and understanding the behavior of composite materials. 
We will build upon existing algorithms 
\cite{Klingenberg1989, Myroshnychenko2005, Cheng2008} to 
efficiently and accurately compute the energy landscape, 
bringing forth the relevant spectral properties and 
critical values of transition phenomena. Small scale experiments, 
involving a few spheres, will isolate the most important 
underlying dynamics in ER fluids as well as provide a 
baseline for numerical simulations. These numerics will 
be extended to include time evolution, providing a full 
scale simulation of the development of chains or fractal 
networks in ER fluids, or the connectedness transition in sea ice.
In conjunction with other projects, we plan to involve the
participants of this study with hands-on experience in
field investigations of sea ice. This  
provides a unique component to this proposal which will
expose the students in particular to the role of
phase transitions in important areas of application,
such as climate change in the polar regions.














\section{Mathematics of transport in two phase media}

\subsection{The analytic continuation method and orthogonal polynomials}
Let us now give a brief overview of the analytic 
continuation method \cite{Bergman:PRL-1285,Milton:APL-300,Golden:CMP-473}. 
The mathematical objects arising here lead
to connections with statistical mechanics and random matrix theory.
We consider a two phase random medium in all of $\mathbb{R}^d$.
Let $\epsilon(\mathbf{x},\omega)$ be the local complex permittivity,
a (spatially) stationary random field
in $\mathbf{x}\in \mathbb{R}^d$ and $\omega \in \Omega$, where $\Omega$ is the
set of all realizations of our random medium. We assume  $\epsilon(\mathbf{x},\omega)$
takes the values $\epsilon_1$ and $\epsilon_2$, and write
$\epsilon(\mathbf{x},\omega) = 
\epsilon_1 \chi_1 (\mathbf{x},\omega) + \epsilon_2 \chi_2(\mathbf{x},\omega)$,
where $\chi_j$ is the characteristic function of medium $j
=1,2$, which equals one for all $\omega \in \Omega$
having medium $j$ at $\mathbf{x}$, and zero otherwise. 



Let $\mathbf{E}(\mathbf{x})$ and $\mathbf{D}(\mathbf{x})$ be the stationary
random electric and displacement fields, related by 
$\mathbf{D}= \epsilon \mathbf{E}$. 
With $\langle \mathbf{E} \rangle=\mathbf{E}_0$, a uniform applied field,
they satisfy \cite{Golden:CMP-473}
\begin{align}                                         \label{gold1_1}  
   \nabla \times \mathbf{E}&=0, \quad
   \nabla \cdot \mathbf{D}=0,\quad \quad
   \mathbf{E}=\mathbf{E}_0+\mathbf{E}_f,    
\end{align}
%
where $\mathbf{E}_f$ is the fluctuating field with mean zero about 
$\mathbf{E}_0$, $\langle \mathbf{E}_f \rangle = 0$, 
and $\langle \cdot \rangle$ 
means ensemble average over  
$\Omega$ or spatial average over all of ${\mathbb{R}}^d$. 
We write 
$\mathbf{E}_0=E_0\mathbf{e}_k$, where
$\mathbf{e}_k$ is a unit vector, for some 
$k = 1, \ldots, d$. 
The effective complex permittivity tensor 
$\beps^*$
is defined as
\begin{equation}                                    \label{gold1_1b}
    \langle \mathbf{D} \rangle=  \beps^* \langle \mathbf{E} \rangle.
\end{equation}
For simplicity, we focus on one diagonal coefficient
$\e^* = \e^*_{kk} =  \langle \epsilon E_k \rangle$.
Central to our studies is the system energy given by 
$\frac{1}{2}\;\langle\mathbf{D}\cdot\mathbf{E}\rangle$ \cite{Jackson}. 
A key variational calculation \cite{Golden:CMP-473}
gives the energy constraint 
$\langle \mathbf{D} \cdot \mathbf{E}_f \rangle=0$,
so that the energy may be written as
%
\begin{equation}                                    \label{gold1_2}
    \frac{1}{2}\;\langle\mathbf{D}\cdot\mathbf{E}\rangle
      =\frac{1}{2}\;\langle\mathbf{D}\rangle\cdot\mathbf{E}_0
        =  \frac{1}{2}\;\mathbf{\epsilon}^*E_0^2.
\end{equation}
%
%For notational simplicity we will focus on one diagonal coefficient
%$\epsilon^* = \epsilon^*_{kk}$, obtained by taking the inner product of 
%$\langle\mathbf{D}\rangle $ with $\mathbf{E}_k$ the $k^{th}$ component of 
%$\mathbf{E}_0$, with $\epsilon^* =  \langle \epsilon \mathbf{E}_k \rangle $,
%and will identify $\mathbf{E}_k$ with $\mathbf{E}_0$.
Due to the homogeneity of $\epsilon^*$, it depends 
only on the ratio $h = \epsilon_1/\epsilon_2$ and we define 
$m(h) =\epsilon^*/\epsilon_2$ \cite{Golden:CMP-473}.

The key step in the method is obtaining
an integral representation for $\epsilon^*$.  
It is more convenient to consider $F(s) = 1 - m(h), \ s = 1/(1-h)$,
which is analytic off $[0,1]$ in the $s$--plane 
\cite{Bergman:PRC-377,Golden:CMP-473}. Then \cite{Golden:CMP-473}
%
\begin{equation}                                     	\label{gold1.eq.3}
  F(s)=1 - \frac{\epsilon^*}{\epsilon_2} =\int_0^1 
   {\frac{d\mu(\lambda)}{s-\lambda}}\;,
   \qquad s = \frac{1}{1-\epsilon_1/\epsilon_2}\;,
\end{equation}
%
where $\mu$ is a positive measure on $[0,1]$. 
This Stieltjes transform of $\mu$ is 
an example of a Herglotz function. 
Applying the spectral theorem \cite{Reed} to the energy 
constraint yields 
\begin{equation}  
\langle E_f^2\rangle=E_0^2\int_0^1 
\frac{\lambda \; d\mu(\lambda)}{(s-\lambda)^2}\;, 
\end{equation}
which in turn leads to an exact, detailed description of all other 
energy contributions within the system 
in terms of Herglotz functions involving $\mu$ \cite{Murphy_Stat_Mech}. The
formula for
$F(s)$ arises from %the resolvent representation  of the electric field, 
$\mathbf{E} = (s + \Gamma \chi_1)^{-1}\mathbf{E}_0$, obtained
via manipulation of (\ref{gold1_1}),
where $\Gamma  = \nabla (- \Delta )^{-1} \nabla \cdot$,
yielding 
$F(s) = \langle \chi_1 [(s + \Gamma \chi_1)^{-1} 
\mathbf{e}_k] \cdot \mathbf{e}_k \rangle$.
%
The operator $\Gamma$ is a
projection onto curl-free fields, based on
convolution with the 
free-space Green's function for $-\Delta$. 
In the Hilbert space $L^2 (\Omega, P)$ with weight $\chi_1$ in
the inner product, $\Gamma \chi_1$ is a bounded self adjoint
operator \cite{Golden:CMP-473}. 
Therefore the underlying self adjoint operator of 
the system is $\chi_1 \Gamma \chi_1$ satisfying 
$\parallel \chi_1 \Gamma \chi_1 \parallel \; \leq 1$.
Formula~(\ref{gold1.eq.3}) 
involves a spectral representation of the resolvent
$(s + \Gamma \chi_1)^{-1}$, where $\mu$ is a spectral measure of the
family of projections of $\Gamma \chi_1$, in the 
$\langle \mathbf{e}_k,\mathbf{e}_k \rangle$ state \cite{Golden:CMP-473,Reed}.
A key feature of \eqref{gold1.eq.3} and \eqref{gold1_2} 
is that the parameter information in $s$ and $E_0$ is {\it separated} from 
the geometry of the composite 
encapsulated in $\mu$, through its moments 
$\mu_i=\int_0^1 \lambda^i d\mu(\lambda)$, which depend on the
correlation functions of the medium \cite{Golden:CMP-473}.
For example, $\mu_0=p_1$, the volume fraction of $\epsilon_1$.

Pad\'e approximants $F_n(s)$ to the Stieltjes 
function $F(s)$ can be written as $F_n(s)=P^{[1]}_n(s)/P_n(s)$, a 
ratio of polynomials with exponential convergence as $n \rightarrow \infty$ 
\cite{Assche}. The denominator is a member of a family of polynomials orthogonal 
with respect to $\mu$ \cite{Ismail}. The recursion relation satisfied by this
family is given by the unitary transformation of $\chi_1\Gamma\chi_1$ to 
tridiagional Jacobi form $J$ \cite{Gautschi}. The roots of the denominator
are eigenvalues of the $n\times n$ principal sub-matrix $J_n$ of $J$ \cite{Gautschi}.
This polynomial family satisfies a minimum principle 
demonstrating that eigenvalue distributions,
corresponding to spectral representations of geometric 
configurations, minimize logarithmic 
potentials. Indeed, the large $n$ limit of the normalized counting measure of eigenvalues, 
$\nu_n=\frac{1}{n}\sum_i \delta_{\lambda_i}$, is the unique equilibrium measure 
of the logarithmic energy supported on $\sigma=supp(\mu)$ (not to be confused with 
conductivity), 
$\mathcal{E}_\sigma[\nu]=\min_{m}\mathcal{E}_\sigma[m]$, 
where $\nu_n \rightarrow \nu$ weakly \cite{Assche,Pastur} and
%
\begin{align}
\mathcal{E}_\sigma[m]=-\int_\sigma U(\lambda;m) dm(\lambda),
\quad U(\lambda;m)=\int_\sigma \log|\lambda-\lambda^\prime|dm(\lambda^\prime),
\end{align}
%  
varied over all probability measures supported on $\sigma$.  The numerator polynomial is also 
a member of a family of orthogonal polynomials \cite{Assche:JCAM-1991-237,Gautschi}. 
%The measure of orthogonality, $d\mu^{[1]}$,
%is not known in general \cite{Assche:JCAM-1991-237,Gautschi} although its Stieltjes transform 
%satisfies $F(s;d\mu)\left[\frac{\mu_2}{\mu_0}F(s;d\mu^{[1]})+\frac{\mu_2}{\mu_0}\right]=F(s;\lambda d\mu)$. 
%The above mentioned Herglotz representations of energy partitions demonstrate 
%this measure is determined by a balance of energy stored in the fields $\chi_1\mathbf{E}$ and 
%$\chi_1\mathbf{E}_d$ making it a central physical 
%object of mathematical importance \cite{Murphy_Stat_Mech}. 




\subsection{Critical behavior in the Ising model and random resistor network}


In order to provide the basis for our statistical mechanics
framework for composites, 
we briefly review the 
nearest neighbor Ising model of a ferromagnet in a field $H$ and 
at temperature $T$ \cite{Thompson-88,Christensen:2005,Golden:JMP-5627}. 
Our treatment will also expose the connection
to the analytic continuation method.
We consider
a finite box $\Lambda \subset \mathbb{Z}^d$ containing $N$ sites.
At each site there is a spin
variable $s_i$ which can take the values $+1$ or $-1$. 
We consider a Hamiltonian with ferromagnetic pair interaction 
$J \geq 0$ between
nearest neighbor pairs $<i,j>$,
\begin{equation}                                                                            
{\cal{H}}_\omega=-J\sum_{<i,j>}{s_is_j}\;-\;H\sum_i{s_i}\; ,                 
\end{equation} 
for any configuration $\omega \in \Omega=\{-1,1\}^N$ of the spins. 
The canonical partition function
$Z_N$ is given by
\begin{equation}                                                        
     Z_N(T,H)=\sum_{\omega \in \Omega}\exp(-\beta {\cal{H}}_\omega) = 
         \exp(-\beta N f_N)  ,      
\end{equation}
where  $\beta=1/kT$, $k$ is Boltzmann's constant, 
and $f_N(T,H)=\frac{-1}{\beta N}\log{Z_N(T,H)}$ 
is the free energy per site.
We are interested in the infinite volume limit 
$f(T,H) = \lim_{N\->\infty} f_N(T,H)$.  
Then the average magnetization                                      
$
M(T,H) =\lim_{N\->\infty}     
\frac{1}{N}\langle \sum_{i=1}^{N}{s_i} \rangle,
$
where $\langle \cdot \rangle$ in this context
denotes averaging over $\omega \in \Omega$ with
Gibbs weights, can be expressed in terms of the
free energy as 
\begin{equation}     \label{mag_free}
M(T,H)=-\frac{\partial f}{\partial H}\;.
\end{equation}
The magnetic susceptibility $\chi(T,H)$ is given by
$\chi(T,H)=\frac{\D M}{\D H}=-\frac{{\D}^2f}{\D H^2}\geq0.$  
When $H=0$, $M(T) \sim(T_c-T)^{\beta}$ as $T\->{T_c}^-$, and 
$\chi \sim (T-T_c)^{-\gamma}$ as $T\->{T_c}^+$. 
Along the critical isotherm $T=T_c$, 
$M(H)\sim H^{1/ {\delta}}$ as $H\->0^+$.
%where $\delta$ and the other exponents have different numerical 
%values  from their analogues in transport or the classical
%percolation model.









The correspondence  between transport in composite media
and statistical mechanics is based on the following 
remarkable property of the partition function.
In 1952, Lee and Yang \cite{Lee:PR-411,Ruelle}
found that the zeros of the partition function of the Ising
ferromagnet (or lattice gas) lie on the unit circle in the $z-$plane, 
where $z=exp(-2 \beta H)$ is the ``activity.''
%$\beta=1/kT$, and $k$ is Boltzmann's constant.
Then %$Z_N$ can be expressed as
\begin{equation}     \label{partition_function}
 Z_N(z) = a_N\prod_{n=1}^N{(z-z_n)}, \;\;\;\; |z_n|=1.
\end{equation}
Equivalently, the Lee-Yang zeros 
lie on the imaginary axis in the $H-$plane. For $T>T_c\;$,
there is a gap $\theta_H$ in these zeros 
around $H=0$ \cite{Ruelle:1971:PRL}, 
which
collapses as $T\->{T_c}^+$, with $\theta_H \sim (T-T_c)^{\Delta}$.
In the $T-$plane
the situation is more complicated,
although in $d=2$ for $H=0$, ``Fisher's zeros'' 
lie on two circles in the complex $v=\tanh (2 \beta J)$ plane
\cite{Fisher-65},  where $J$ is the interaction strength.










The Lee-Yang Theorem leads to an elegant analytic structure for 
the Gibbs free energy per site $f(T,H)$. For finite volume,
%$f_N=-\frac{1}{\beta N}\log{Z_N}$, or
%\begin{equation}        \label{free_energy_N}
$f_N = -\frac{1}{\beta} \sum_{n=1}^N \frac{1}{N} \log{(z-z_n)} -
       \frac{1}{\beta N} \log{a_N}.$
%\end{equation}
In the thermodynamic limit as $N \rightarrow \infty$, 
\begin{equation}        \label{free_energy}
 f(T,H) = -\frac{1}{\beta} \int_0^{2\pi} \log{(z-e^{i\theta})}d\nu(\theta) -
                      2d\beta J,
\end{equation}
where $\nu$ is a measure on the circle describing the limiting distribution
of Lee-Yang zeros, and $d$ is dimension. 
The integral is one half the electrostatic potential of
the corresponding distribution of point charges on the unit circle.  
The minimization of closely related sums and integrals on manifolds,
and their relation to many problems in mathematics and applications is 
surveyed in \cite{Hardin:NAMS-2004}.
Baker \cite{Baker:PRL-990} used this logarithmic potential representation
to show that $M$ has the following
analytic structure in $\tau = \tanh(\beta H)$,
\begin{equation}              \label{baker_magnetization}
  M(\tau) = \tau + \tau(1-\tau^2)G(\tau^2), \;\;\; \qquad 
   G({\t}^2)=\int_0^{\infty}\frac{d\psi(y)}{1+{\t}^2 y}\;,
\end{equation}
where $G$ is a Stieltjes (or Herglotz) function of $\tau^2$,
and $\psi$ is a positive measure which for $T>T_c$ is supported 
only in $[0,S(T)]$, where $S(T) \sim (T-~T_c)^{-2\Delta}, T\->{T_c}^+.$
This integral representation was used to obtain the 
classical (two-parameter) scaling relations of statistical  mechanics,
$\beta = \Delta - \gamma$ and $\delta = \Delta / (\Delta-\gamma)$ 
\cite{Gaunt:PRB-1184,Baker-1990}, 
%(which are satisfied by
%the mean field exponents 
%$\beta=1/2, \gamma = 1, \delta = 3, \Delta = 3/2$, 
%and the exact exponents for $d=2, 
%\beta=1/8,  \gamma = 1 \frac{3}{4},  \delta = 15, 
%\Delta=1 \frac{7}{8}$),
as well as {\it Baker's inequalities}
$\gamma_{n+1} - 2\gamma_n + \gamma_{n-1} \geq 0$ for the critical
exponents $\gamma_n$ of the higher field derivatives of the free energy $f$,
or equivalently, of the moments $\psi_n$ of $\psi$,
$\psi_n\sim(T-T_c)^{-\gamma_n}, T\->{T_c}^+$
with $\gamma_0 = \gamma$ \cite{Baker:PRL-990}. 
%The sequence $\gamma_n$
%is actually linear in $n$, $\gamma_n=\gamma+2n\Delta , n\geq0$, 
%with constant gap 
%$\gamma_i-\gamma_{i-1}=2\Delta$ \cite{Fisher:RPP-615,Baker:PRL-990}.



It was observed in \cite{Golden:PRL-3935} that the
Lee-Yang-Baker critical theory outlined above applies to classical transport
problems such as electrical conductivitiy and
permittivity, and magnetic permeability.
In particular, we obtained a new Stieltjes
integral representation 
\begin{equation}
  m(h) = 1 + (h-1) g(h), \qquad \;\;\; 
  g(h)=\int_0^{\infty}\frac{d\eta(y)}{1+h y}\;,
\end{equation}
where $\eta$ is a positive measure depending only on the geometry
of the composite and $y=\lambda/(1-\lambda)$, which is the
direct analogue of Baker's formula (\ref{baker_magnetization}), 
and is a different form of 
formula (\ref{gold1.eq.3}).
For the percolation models such as the
random resistor network, we assume that 
for $p>p_c$, there is a gap 
${\theta_h}\sim(p-p_c)^\Delta$ in the support of $\mu$ around 
$h=0$ or $s=1$, (or that any spectrum in 
this region does not affect power law behavior, which is discussed
in detail in \cite{Golden:PRL-3935} and is supported by 
\cite{Clerc:AP-191,Day:JPCM-96,Bruno:PRSLA-353}), 
so that $\eta$ is supported 
only in $[0,S(p)]$, where $S(p) \sim (p-~p_c)^{-\Delta}, \; p\->{p_c}^+.$
Then, the methods developed by Baker 
\cite{Baker:PRL-990,Gaunt:PRB-1184,Baker-1990}
to analyze the critical behavior of the Ising model can be
modified to obtain scaling relations for the
critical exponents for transport,
$ t = \Delta - \gamma,\;
\delta =\Delta/(\Delta-\gamma).$
The relations hold for lattice and continuum models, even though
continuum models can exhibit non-universal behavior
\cite{Halperin:PRL-2391,Stauffer-92,Torquato:RHM-02}.
The critical
exponents $\gamma_n$ of the higher derivatives of $m(h)$ at $h=0$,
or the moments $\eta_n$ of $\eta$,
with $\eta_n(p) \sim (p-p_c)^{-\gamma_n}$,
also satisfy Baker's inequalities.
%Furthermore, like the Ising model,
%$\gamma_n$ for transport is linear in $n$, $\gamma_n=\gamma+(n-1)\Delta , n\geq1$, 
%with constant gap $\gamma_i-\gamma_{i-1}=\Delta$, consistent with
%the absence of multifractal behavior for the bulk conductivity \cite{Stauffer-92}. 




For a finite sample of a resistor network (or continuum model)
with $N$ resonances $\lambda_n\in[0,1]$, or poles of $F(s)$,
we can construct an analog of (\ref{partition_function}) 
for transport,
%\begin{equation}
${\Z}_N(s) =\prod_{n=1}^N{(s-\lambda_n)}$, 
%\end{equation}
whose zeros are distributed according to a
normalized counting measure by
$\nu_N=\frac{1}{N}\sum_{n=1}^N \delta_{\lambda_n}$.
Then we define a
finite volume ``Gibbs free energy'' \cite{Golden:JMP-5627} via
$
\Phi_N(s)=-\frac{1}{N}\log{\Z}_N(s)=-\int_0^1 \log(s-\lambda)d\nu_N(\lambda).
$
In the infinite volume limit \cite{Assche},  
\begin{equation}        \label{free_energy}
\Phi(s)=-\int_0^1 {\log(s-\lambda)d\nu(\lambda)}.
\end{equation}
It is a standard result of logarithmic potential theory that $\nu$ is the unique equlibrium measure of the ``electrostatic energy'' given by $\mathcal{E}[\nu]=\inf_m\mathcal{E}[m]=\inf_m\left( -\int_{\sigma\times\sigma} \log|\lambda^\prime-\lambda|dm(\lambda)dm(\lambda^\prime) \right)$ varied over all probability measures $m$ supported on $\sigma=supp(\nu)$ 
\cite{Pastur,Assche}.
We expect this minimum principle to aid in the characterization
of ground states and equilibrium configurations of composite systems.


%It is a standard result of logarithmic potential theory that $\Phi(s)=\inf_m\left[ -\int_\sigma \log|s-\lambda|dm(\lambda) \right]$ varied over all probability measures $m$ supported on $\sigma=supp(\nu)$ \cite{Pastur,Assche}.

%******Logarithmic capacity - Ben********







In \cite{Golden:JMP-5627,Golden:PRL-3935} 
the zeros of the transport partition function ${\Z}_N(s,p)$ in the complex
$p-$plane were considered. They are the analogs of Fisher's zeros 
\cite{Fisher-65} in the complex temperature plane.
For the $d=2$ random resistor network we found that the
percolation threshold $p_c=\frac{1}{2}$ can be characterized as 
an accumulation point of
its zeros.  
%It was also proven \cite{Golden:JMP-5627} that 
%for {\it infinitely interchangeable media}, (see, e.g. \cite{Bruno:JSP-365}), 
%which include lattice and cell materials in the 
%continuum, such as the random checkerboard,
Also, $m(p,h)$ is analytic in $p$ for $|h-1|<1$ in a 
domain which includes [0,1],
rigorously establishing the existence of a gap $\theta_p$ in the zeros 
in the $p-$plane around $p_c$. With $\theta_p \sim (p-p_c)^{\Delta_p}$,
${\Delta_p} \leq 1$.










\subsection{Numerical Modeling of Transport in Two Phase Composites}
%
Numerically simulating particles in a host medium accurately and efficiently 
has been a subject of longstanding as well as recent interest in the applied
mathematical community. Computationally, the spectral properties of transport 
phenomena can be investigated using various techniques, 
such as molecular dynamics-like formulations, 
fast multipole methods, and finite element methods.

When an electric field is applied to an electrorheological fluid, 
chain or fractal network structures result from the induced dipoles 
on the particles.
Dynamically, the formation of chains has been obtained by considering 
the forces on each particle due to the electric field, 
induced dipoles of other particles, and the surrounding fluid  
\cite{Klingenberg1989}. The equation of motion for each particle, 
when integrated in time, allows aggregates to form into a percolating
cluster. Chain growth has also been modeled using 3D Monte 
Carlo simulations \cite{Weis1993}, and by approximating the chain 
as a prolate spheroid \cite{Wang2001}. Beginning with an initial 
random configuration in a continuum, the development of fractal 
patterns was numerically simulated by dynamically calculating the 
effects of particle aggregation on the applied electric field 
\cite{Kun1998}. This model, while showing agreement with the 
fractal dimensions obtained experimentally, is restricted to 
low particle concentrations.

More recently, finite element methods have been used to accurately 
calculate the effective permittivity of random inhomogeneous media. 
Configurations of dielectric inclusions in a host matrix have been 
considered in the 2-d case with circular, non-overlapping particles 
\cite{Myroshnychenko2005}. In a different study, random placements 
of inclusions in a 3-d block of cubic cells allowed a wide variety 
of interconnecting inclusion geometries \cite{Cheng2008}. 
For both of these models, the effects of the anisotropy in a 
single configuration was eliminated by ensemble averaging. 
The energy $W$ in the parallel plate capacitor of area $S$ and 
separation distance $h$ is related to the effective permittivity 
$\epsilon^*$ by the %free space permittivity $\epsilon_0$ and 
change $u_0$ in the applied potential between the plates,
%\begin{equation}
$W=\frac{1}{2}\epsilon^* \frac{S}{h}u_0^2$.
%\end{equation}
This energy can also be computed as an integral over 
the volume $V$ of the capacitor,
\begin{equation}W=\frac{1}{2} 
\int_{V} \epsilon(\mathbf{x}) 
[\nabla u(\mathbf{x})]^2 d\mathbf{x}.
\end{equation}
%which takes into account the local dielectric constant 
%$\epsilon(\mathbf{x})$ and local electrostatic potential 
%$\phi(\mathbf{x})$. 
The electrostatic potential $u(\mathbf{x})$ is
found by solving the boundary value problem, 
$\nabla \cdot [\epsilon(\mathbf{x}) \nabla u(\mathbf{x})]=0$, 
where $u=0$ on the bottom plate and 
$u=u_0$ on the top plate. With this formulation, 
along with commercial finite element solvers, we can compute 
the spectral properties of specified, static configurations 
of ER fluids or sea ice structures.
We can then build a comprehensive, accurate, and efficient 
numerical tool to compute the dynamic or static energy and 
effective permittivity of two phase composites.







\section{Statistical Mechanics of Homogenization for Composites}




We now focus on the {\bf proposed research}, the questions
we plan to address, and the tools we will develop to carry out
our investigations. 
There has been a persistent and fundamental puzzle in 
electrorheology. We expect that
a statistical mechanical understanding
of the effective properties of composites will
help resolve it, and may yield 
an entirely new way of looking at field-induced 
percolation transitions.

Consider a system of identical millimeter-sized 
dielectric spheres placed on a flat, rectangular 
surface bounded on two sides by parallel electrodes.  
The friction between the surface and the spheres 
can be made very small so that its effect can be 
ignored.  When a small voltage is applied across 
the electrodes, the dielectric spheres will move 
and form necklaces roughly aligned to 
the electric field.  But percolation, in the sense 
that touching spheres form a connected network 
bridging the two electrodes, will occur only at 
a sufficiently large voltage that is approximately 
repeatable if the system is large enough.    
If instead we have a 3-d suspension 
of identical, ten micron-sized dielectric spheres 
in an insulating fluid, the same phenomenon occurs.
The large size of the spheres means that the 
temperature effect can be ignored. That is, the field required
for percolation is observed to be far larger than would be required
to overcome thermal fluctuations, since
for a typical field on the order of 1,000 V/mm,
$\alpha r^3 E_0^2/kT \gg 1 $, where $r$ is the sphere 
radius and $\alpha=(\epsilon_1 - 1)/(\epsilon_1 + 2)$ is the 
Claussius-Mossotti factor for spheres of dielectric constant 
$\epsilon_1$.  

The puzzle is: what determines the electric field 
at which the percolation transition occurs?  
Posed in another way: where does the threshold field come from?
If one looks at the problem from the energy point of view, 
there seems to be no energy scale in the problem apart 
from the electric field, if $kT$ can be discounted.  
So what determines the threshold field value?  
Theoretical calculations show that under any electric 
field, the body-centered-tetragonal (BCT) arrangement 
of the spheres constitutes the lowest energy 
microgeometry.  Hence the spheres should move into such an 
arrangement under any value of the applied 
electric field.
However, one observes that a threshold 
value for the field has to be exceeded before such 
chains/columns can span the two electrodes, and
long-range order is achieved.

A potential solution 
lies in the observation that there could be many 
metastable configurations for the spheres under 
a moderate electric field, and the percolation 
transition point occurs when the variance of 
the energy landscape vanishes, or becomes on 
the order of $kT$. Alternatively, it may be that 
a certain field strength is required to overcome
the ``typical'' energy barrier, which may itself 
evolve with increasing field.
The observation that there 
could be many metastable configurations is 
supported by simulations 
showing that the movement of the spheres would 
stop at intermediate configurations before reaching 
percolation. Stopping implies a force balance, 
or a metastable configuration.  
But metastability also implies barriers in the energy landscape.  
Therefore, analysis of the energy terrain,
its variance, and the barriers associated with metastable
states should yield fundamental insights into the problem.
Such an analysis will open new perspectives 
on the ER fluid, and may serve as a paradigm 
for other geometric transitions. 










To gain initial insight into this type of behavior,
we developed a 2-d square lattice dipole model of a dielectric 
ER fluid, with its evolution based on Monte Carlo simulated annealing.
The energy is calculated from
a simplified Hamiltonian incorporating
interactions of the dipoles with the applied field and their neighbors.
The probability that a given configuration evolves to one 
where a particle moves is determined by the associated Gibbs factor. 
The model captures the qualitative features of 
low field, disordered states and high field, columnar 
states \cite{Murphy_Stat_Mech}, whose evolution  
as the field is increased
is shown in Figure \ref{simulated_annealing}.
We find in our simulations the very interesting phenomenon
discussed above.
Graphs of the system energy show an overall descent
as the field is increased, but with small hills that 
must be scaled before descending further. Comparison with
configurations indicates that the system
temporarily gets stuck in slightly higher energy configurations,
until the field is strong enough to ``kick'' it into
an overall lower energy state. We believe this type of behavior
is at the heart of the existence of a critical field.
We have seen a similar phenomenon - a hilly energy landscape - 
in another type of simulation where
we let the system evolve to a new state if it lowered 
the logarithmic integral in (\ref{free_energy}).
In these simulations with a relatively small number
of particles we observed chain formation as the field increased
beyond a critical point. 
Through more detailed numerical simulations
we will extend these models 
to continuous systems of disks and
spheres, and 
obtain the energy of a configuration much more precisely by direct 
calculation of the electric field,
providing a more accurate description of ER fluid 
evolutions and the energy landscape.



\begin{figure}                \centering
\includegraphics{spectral_measures_ben.eps}
   \caption{A simulated annealing model of 
microstructural evolution in electrorheologial fluids.
Long range order develops as the applied electric field is 
increased (left to right) and exceeds a critical value $-$ 
similar to sea ice as temperature increases. 
The system tends to metastable states determined 
by the Gibbs factor associated with a system Hamiltonian.  
As the field strength increases, 
these metastable states destabilize and 
evolve toward columnar, checker board ground states.}
      



    \label{simulated_annealing}
        
\end{figure}









When an electric field is applied to a dielectric ER fluid, 
the particles become polarized creating interactive forces 
which change the underlying geometry. As the geometry 
evolves, the associated spectral measure evolves 
accordingly. 
Generically, the equilibrium 
distribution of a statistical mechanics model is given 
by the global minimum of the associated Gibbs-Helmholtz 
free energy, and during phase transition, the 
analytic properties of this free energy change 
suddenly \cite{Robertson,Baker-1990,Ruelle}. In the ER fluid
this picture may be complicated by temporary, local minima.
One of our main goals is to characterize 
equilibrium configurations and phase 
transitions in ER fluids and other composite 
systems by the minima and analytic 
properties of the free energy, and to obtain powerful
representations which facilitate this analysis.






Toward this goal, we present here a statistical mechanics 
model of two phase composites, and
ER fluids in particular.
%which can be extended to the 
%case of metallic spheres by letting $s\to0$.
We consider a finite box $B_V \subset \mathbb{R}^d$ of 
volume $V$ containing a number density $n$ of spheres 
in an electric field $\mathbf{E}_0$ and at temperature 
$T$. One obtains the Hamiltonian from the system 
energy by exclusion of the contribution associated with 
the external source, and including energy
stored within the system and which performs work on 
the system. Herglotz representations facilitate 
an especially simple energetic analysis arriving at the 
following system Hamiltonian \cite{Murphy_Stat_Mech,Robertson}
%
\begin{align} 													\label{eq:ER_fluid_Hamiltonian}
\mathcal{H}_{\omega}(E_0,s)=
- \frac{1}{2}P_\omega^*(E_0,s) E_0=- \frac{1}{2n}
\epsilon_2 E_0^2\chi_{\omega}^*(s)
\end{align}
%
for any configuration $\omega\in\Omega$. Here 
$\epsilon_2$ is the permittivity of the liquid 
host, $P_\omega^*$ plays the role of the polarization, and 
$\chi_{\omega}^*(s)=F_{\omega}(s)-p_1/s$ plays the 
role of the electric susceptibility with an integral 
representation 
\eqref{gold1.eq.3} with $\omega$ dependent measure 
\cite{Murphy_Stat_Mech}. The canonical partition 
function $Z_V$ is given by
%
\begin{align}
	Z_V(T,E_0,s)=\sum_\omega 
\exp(-\beta \mathcal{H}_{\omega}(E_0,s))=\exp(-\beta V f_V),
\end{align}
%
where %$\beta=1/kT$ and $k$ is Boltzmann's constant, 
$f_V(T,E_0,s)=\frac{-1}{\beta V}\log{Z_V(T,E_0,s)}$ 
is the free energy per unit volume. We are interested 
in the infinite volume limit 
$f(T,E_0,s)=\lim_{V\to \infty}f_V(T,E_0,s)$ such 
that the particle density $n$ remains constant. 
The average polarization, with $\langle \cdot \rangle$ 
denoting averaging with respect to Gibbs factors, 
is given in terms of $f$ by 
%
\begin{align}
P^*(T,E_0,s)=-\frac{\partial f}{\partial E_0} =
\frac{\epsilon_2 E_0}{n} \left(\langle F_\omega(s)\rangle - \frac{p_1}{s}\right)=
\frac{\epsilon_2 E_0}{n} \left(F(T,E_0,s)-\frac{p_1}{s}\right).
\end{align}
%
The average electric susceptibility is given by 
$-\partial^2 f/\partial E_0^2=
\frac{\epsilon_2}{n}\left(F(T,E_0,s)-\frac{p_1}{s}\right)
- \frac{\beta}{n}(\epsilon_2 E_0)^2\text{Var}(F_\omega)$
\cite{Murphy_Stat_Mech,Chandler}, which determines the 
local convexity properties of $f$ in $E_0$. The possible
competition in sign here may be related to the issue of 
metastability. 
%as the $p_1/s$ term vanishes in the limit. 
%The average electric susceptibility is given by 
%$-\partial^2 f/\partial E_0^2=
%\partial P^*/\partial E_0=
%(\epsilon_2 F(T,E_0,s)-\beta 
%(\epsilon_2 E_0)^2\text{Var}(F_\omega))/n$ 
%\cite{Murphy_Stat_Mech} which determines the 
%convexity properties of $f$ in $E_0$. 
This is a central object in the characterization of $E_c$. 
Implications of the fluctuation dissipation theorem 
\cite{Chandler} regarding the average susceptibility 
are key in this direction of research. 

The polarization $P^*(T,E_0,s)$ is called a {\it response 
function} since it determines the response of the system 
to changes in system parameters. Other response 
functions are related via
\cite{Murphy_Stat_Mech} 
%
\begin{align}
df=-SdT-P^*dE_0-\Psi ds, 
\end{align}
%
where $S$ and $\Psi$ denote the entropy and 
contrast response functions, respectively.

From the first law one can relate changes in response 
functions due to variations in system parameters. 
These so-called Maxwell's relations are obtained by 
equating mixed derivatives of $f$
\cite{Robertson,Baker-1990}, e.g., 
%%
\begin{align}
	-\frac{\partial^2 f}{\partial s E_0}=
\frac{\partial P^*}{\partial s}=
\frac{\partial\Psi}{\partial E_0}. 
\end{align}
%%
A surprising consequence of Maxwell's relations 
is that the covariance of the Hamiltonian and the 
contrast derivative is identically zero independent 
of the details regarding the evolution to equilibrium, 
geometry class, and parameter values. 
This universal constraint on the system implies 
%%
\begin{align}
\text{Cov}\left( F_{\omega}, F_{\omega}^\prime\right)
        =\langle   \left( F_\omega  -  
       \langle  F_\omega  \rangle        \right)
         \left( F_{\omega}^\prime -  \langle  F_{\omega}^\prime  
           \rangle   \right) \rangle \equiv 0
\end{align}
%
\cite{Murphy_Stat_Mech}, where prime denotes differentiation 
with respect to $s$. 
This constraining equation dramatically simplifies 
the analysis in the 
$s$ variable, yields explicit analytical 
representations of the response functions $S$ and 
$\Psi$, and connects $\Psi$ with the energy 
contained in $\chi_1 E^2$ 
\cite{Murphy_Stat_Mech}. Maxwell's relations 
also demonstrate that fluctuations in macroscopic 
parameters are governed by a logarithmic functional 
$\Phi(T,E_0,s)\equiv\ln{F(T,E_0,s)}$.
%, e.g. 
% $ -\partial s/\partial T=\partial S/\partial\Psi=(\partial\Phi/\partial T)(\partial\Phi/\partial E_0)^{-1}$.
%%\end{align}
%%
%These relations constitute a
%detailed balance of energy in the system, which can be shown by perturbation analysis \cite{Murphy_Stat_Mech}.
This logarithmic functional takes on a more suggestive 
form when we consider Pad\'{e} approximants of $F(T,E_0,s)$,
$F_m= P_m^{[1]}/P_m$, yielding 
$\log{F_m}=\sum_i  (\log(s-\lambda_i^{[1]}) -  \log(s-\lambda_i))$,
where $\lambda_i^{[1]}$ and $\lambda_i$ are 
roots of the polynomials $P_m^{[1]}$ and $P_m$, respectively.
Denoting the infinite volume root densities by 
$\rho^{[1]}$ and $\rho$ \cite{Assche,Ismail} 
%
\begin{align} \label{eq:Log_Potential}
  \Phi(T,E_0,s)=-\int_0^1\ln(s-\lambda)\langle 
d\rho\left( \lambda \right)-
d\rho^{[1]}\left( \lambda \right)\rangle
=U\left( s;\langle\rho\rangle \right)-
U(s;\langle\rho^{[1]}\rangle ),
\end{align}
%
in the infinite volume limit. 
Thus, well-developed techniques of logarithmic 
potential theory \cite{Staff_Totik,Deift} can be used to study the energetic
processes described by $\Phi$ in \eqref{eq:Log_Potential}. 
A key question is whether or not a minimum principle 
is associated with $\Phi$.
Our numerical simulations of a discrete approximation 
of the 2--d composite checkerboard have demonstrated 
that minimization of the logarithmic 
potential $U(s;\mu)$ yields chain-like ground 
states for dielectric inhomogeneities, and fractal 
net ground states for metal inhomogeneities.
This minimum principle is also of key interest 
in the characterization of ground states in ER fluids.

The above framework mathematically 
parallels the Ising model, providing guiding 
intuition for our future investigation of ER fluids. 
For example, the onset of critical behavior in the 
Ising model is caused by the collapse of a gap in the
partition function zeros as 
$T\to T_c$ . Likewise, non-analytic behavior of 
the free energy in binary composites is caused 
by the collapse of spectral gaps in $F_\omega(s)$ 
as $p \rightarrow p_c$ for high contrast materials, 
leading to non-analytic behavior in 
$F_\omega(s)$ \cite{Golden:PRL-3935}. 
Furthermore, just as the separation of parameters 
in the Ising model gives a polynomial representation 
of the partition function in 
$\exp(-2\beta H)$, the parameter separation in binary 
composites gives a fractional polynomial representation 
of the partition function in the variable 
$\exp(-\beta \epsilon_2 E_0^2/2n)$. One of the 
fundamental questions which must be addressed is 
whether or not this representation yields a 
factorization of the partition function which 
leads to a logarithmic integral representation 
of the associated free energy, as in the Ising model. 
Recent work suggests the family of polynomials 
orthogonal to $\mu$, complete in the Hilbert space 
$L^2(\mu)$, may provide a natural basis to expand 
the partition function in order to obtain a 
logarithmic integral representation of the free 
energy. If such a representation exists, 
well-developed methods of logarithmic 
potential theory \cite{Staff_Totik,Deift} could 
also be used to characterize ground states of the 
ER fluid as minima of the potential.












\section{Random Matrix Theory and Phase Transitions in Composites}





\begin{figure}                \centering
\includegraphics{Figure_Unfolded_Spacing_Density_SF.eps}
   \caption{Spectral statistics of the 2-d bond RRN. (a) The nearest neighbor eigenvalue spacing distribution for various volume fractions $p$, unfolded to have unit local mean. The insulating regime ($p<p_c=0.5)$ is characterized by Poisson-like statistics which transition, with $p$ increasing to 0.5, toward Wigner Dyson (WD)-like statistics in the conducting regime.  (b) The corresponding evolution of the spectral function $\mu(\lambda)$. As $p\to p_c$ gaps in the spectral function collapse leading to critical behavior in $\epsilon^*$.} 

    \label{spectral_measures_ben}
        
\end{figure}






We have found that critical behavior of binary composites 
can be studied using a broad range of mathematical techniques. 
Here we show that random matrix theory (RMT) also provides 
invaluable insight into this problem. A central object 
in homogenization theory of composites is the self adjoint 
random operator 
$\chi_1\Gamma\chi_1$, where 
$\Gamma=\nabla(-\Delta)^{-1}\nabla\cdot$ and $\chi_1$ 
is the characteristic function of type one medium 
\cite{Golden:CMP-467}. 
The action of this random integro-differential operator 
on a bond random resistor network (RRN) in $\mathbb{Z}^d$ is
given by that of a real sparse symmetric random matrix 
\cite{Golden:CMP-467}. The differential operators $\nabla$ 
and $\nabla\cdot$ are given by forward and backward 
difference operators, respectively, 
$(-\Delta)^{-1}$ is the inverse of the lattice Laplacian, 
and $\chi_1$ is a random diagonal matrix of ones and zeros 
\cite{Golden:CMP-467,Murphy_RMT}.  This structure 
allows us to incorporate powerful analytical and 
variational techniques of RMT into our studies of phase 
transitions in disordered binary composites.
 
Techniques of RMT have been used to accurately describe 
phase transitions in many disordered systems,  
including random graphs, quantum chaos, and disordered 
quantum and mesoscopic conductors 
\cite{Efetov,Canali,Muttalib_etal_qRME,Kuhn}. 
The random matrix description of all these systems 
have common universal features. Most well known is 
the transition in eigenvalue statistics from highly 
correlated Wigner-Dyson (WD) statistics toward uncorrelated 
Poisson statistics, as a function of disorder 
\cite{Muttalib_etal,Canali}. The strong eigenvalue 
correlations in WD statistics gives rise to the 
phenomenon of eigenvalue repulsion 
\cite{Mehta:2004:RM,Deift}. 

In order to observe the universal eigenvalue fluctuations 
of different matrices, random matrix theory requires 
spectral unfolding to remove the configuration specific 
local mean spacing. Our numerical simulations 
of the binary bond RRN on $\mathbb{Z}^2$ demonstrate 
\cite{Murphy_RMT} that the eigenvalue 
statistics exhibit this transition as a function of 
volume fraction, as shown in Figure \ref{spectral_measures_ben}. 
Similar behavior was observed in an 
approximately symmetric, random matrix arising from a 
Green's function formalism \cite{Gu:PRE:2002}.
These results illustrate that macroscopic 
electric systems exhibit the same {\it universal} fluctuations 
as mesoscopic electric systems of size $10^{-6}$ m and 
quantum systems nine orders of magnitude smaller 
\cite{Muttalib_etal,Evangelou:PRL:1992}. The reason 
that the statistics cannot complete the transition to 
the WD regime as $p \to p_c$ is that WD statistics are 
independent of the system parameters 
\cite{Muttalib_etal,Canali}, which is not true 
for RRN's near percolation. 

Our simulations 
show that gaps in the eigenvalue distribution collapse 
before those of the spectral function, setting the 
stage for critical behavior. The spectral function 
gaps subsequently collapse, with increasing $p$, giving 
rise to critical behavior at the percolation threshold 
$p_c$. Intuitively, the behavior is symmetric about 
$p=0.5$. For example, the figures for $p$ and $1-p$ 
are identical, up to a scaling factor, when the 
buildup at the spectral endpoints is discarded 
\cite{Murphy_RMT}. For this reason, our figures 
for $p>p_c$ differ from those in \cite{Day:JPCM-96} 
which were found by direct calculation of $\sigma^*$ 
from application of the Stieltjes--Perron inversion 
theorem yielding $\mu$ \cite{Gautschi}. 
We will generalize our software to higher dimensions 
and calculate other important eigenvalue statistics 
\cite{Mehta:2004:RM} in order to gain further insight 
into the implications of this transition in eigenvalue 
statistics and associated change in the spectral 
function. We will also use analytical and variational 
techniques 
\cite{Deift,Muttalib_etal,Evangelou:PRL:1992,Mehta:2004:RM}, 
discussed below, to characterize phase transitions and 
ground-states in binary composite media. The RMT 
of mesoscopic and quantum conductors will serve as an intuitive 
guide in the characterization of phase transitions in 
macroscopic conductors. 






\begin{figure}                \centering
\includegraphics{DLA_Spectral_Function.eps}
   \caption{Evolution of the spectral function for a growing 
   DLA cluster. The plots exhibit a decreasing spectral gap as the cluster percolates to the boundaries. }



    \label{DLA_Spectral_Function}
        
\end{figure}









The probability distribution of symmetric random matrix 
ensembles (RME's) can be expressed solely in terms of 
its eigenvalues $\{\lambda_i\}$
\cite{Mehta:2004:RM,Deift}. Surprisingly, this 
distribution is in the exact form 
\cite{Mehta:2004:RM,Deift,Balian,Dyson} of 
the statistical mechanics ensemble of 
positive charges on a line in two dimensions, in an external confining 
potential $V>0$, competing with logarithmic 
repulsion: $P=Z_N^{-1}\exp(-\beta \mathcal{H}_N)$, where     
%
\begin{align}
  \mathcal{H}_N\{\lambda_1,\ldots,\lambda_{N}\}=
\beta^{-1}\sum_{k=1}^{N}V(\lambda_k)-
\sum_{k>l=1}^{N}\ln|\lambda_k-\lambda_l|.
\end{align}
%
Here $N$ is the size of the matrix, $Z_N$ is 
the partition function, and $\beta$ depends only 
on the symmetry class taking values 1 and 2 for 
real and Hermitian symmetric matrices, respectively
(not to be confused with $\beta = 1/kT$). 
When $V(\lambda)\geq\lambda$ universal WD statistics 
arise, in the infinite volume limit $N\to\infty$, 
which are independent of all system parameters 
\cite{Muttalib_etal_qRME,Canali}. At the onset 
of criticality $V$ violates this inequality 
becoming comparable to $\ln^2{\lambda}$ and 
the statistics deviate toward Poisson 
statistics \cite{Muttalib_etal_qRME,Canali}. 
In the Hermitian ($\beta=2$) case, it has been 
analytically proven that, in this regime, the 
statistics are dependent only on a single 
disorder parameter \cite{Muttalib_etal_qRME}. 
It is expected that the orthogonal ($\beta=1$) 
case also has this feature
\cite{Muttalib_etal}.

In total agreement with statistical mechanics 
\cite{Pastur}, for $V$ satisfying 
$\lim_{|\lambda|\to\infty} V(\lambda)/
\log(\lambda^2+1)=\infty$, the ground state 
distribution of eigenvalues is given by the 
minimum ``electrostatic'' energy. 
This distribution is determined by the 
unique equilibrium measure $\rho$ of the 
associated infinite volume free energy: 
$f[\rho]=\inf_m f[m]=\lim_{N\to\infty}N^{-2}\ln{Z_N}$, 
where 
%
\begin{align}
  \label{eq:GroundState_Measure_min}
  f[m]=
     \int_0^1\left(V(\lambda)-\frac{\beta}{2} 
     U(\lambda;m)\right) dm(\lambda), 
\quad U(\lambda;m)=\int_0^1\ln|\lambda-
 \lambda^\prime|dm(\lambda^\prime).                  
\end{align}
%
Demonstrating, once again, the ground-states of 
a binary composite are characterized by the 
minimum of a free energy. We stress that the 
corresponding variational problem determines both 
the (compact) support $\sigma$ of the measure and 
the form of the measure \cite{Pastur}. Moreover, 
the corresponding Euler-Lagrange equations show 
that, in the infinite volume limit, the confining 
potential $V$ satisfies the singular integral equation 
$V^\prime(\lambda)=\int_{\sigma} d\rho(z)/(\lambda-z)$ 
for $\lambda\in\sigma$ \cite{Pastur}. These key 
properties will be aggressively pursued in our 
studies of phase transitions of binary composites. 

The ground-state distribution of eigenvalues depends 
on the form of the confining potential $V$ although 
fluctuations about this distribution are universal 
\cite{Muttalib_etal_qRME,Canali}. This fact has been 
exploited to obtain extremely detailed information 
about phase transitions in disordered systems using 
powerful analytical methods  
\cite{Muttalib_etal_qRME,Canali,Efetov,Evangelou:PRL:1992,Kuhn}. 
For Hermitian RME's ($\beta=2$) Mehta's orthogonal polynomial 
method has led to exact closed form calculations of all 
relevant statistics and all $n$-point correlation functions 
\cite{Mehta:2004:RM,Muttalib_etal_qRME}. 
This gives an exact, detailed description of phase 
transitions in mesoscopic/quantum conductors described 
by Hermitian random matrices 
\cite{Muttalib_etal_qRME,Canali}. Similar results are 
expected to follow in our orthogonal RME ($\beta=1$) 
\cite{Muttalib_etal,Canali}. In the absence of a 
generalization of Mehta's method to orthogonal RMEs, 
we will explore the well-developed super-symmetry  
\cite{Efetov,Canali} and replica methods   
\cite{Efetov,Canali,Evangelou:PRL:1992,Kuhn} in describing 
phase transitions in eigenvalue statistics of real 
symmetric, sparse random matrices.  There is a natural 
question regarding Pad\'e approximants of 
$F_N(s)= P_N^{[1]}(s)/P_N(s)$. If the numerator polynomial 
zeros, like the denominator zeros, also follow these 
universal eigenvalue statistics, this would imply 
$\mu$ itself has universal features which would 
completely couple physical geometric fluctuations 
to these universal fluctuations. 

In Figure \ref{DLA_Spectral_Function} we show what we believe to be the 
first calculations of the spectral properties of diffusion limited
aggregation (DLA) clusters. In this purely probabilistic 
process \cite{Meakin:1998,Bunde:FDS-91}, a particle
is placed at the origin, and a random walker stops if it lands
on a nearest neighbor site, creating a new cluster with two sites.
A new walker enters and stops if it lands on a nearest neighbor of that
cluster, and so on. The growth of the cluster is quite
similar to the growth of the fractal networks in metal ER fluids.
We plan to explore the similarities between these processes
through RMT calculations, and
the graphs in Figure \ref{DLA_Spectral_Function} represent preliminary
work in this direction. We suspect that there may be
formulations of the two problems, involving the Laplacian  
\cite{Meakin:1998,Bunde:FDS-91},
which will eventually display deep mathematical parallels
between them. In conjunction with our characterizations of
ER fluid ground states as minimizers of logarithmic potentials,
we will explore the characterization of DLA clusters as 
solutions of variational principles, opening a new chapter
in the analysis of this imporant process.







\section{Experiments and Numerical Simulations}
%

Our close colleagues Ping Sheng and Weijia Wen at the 
Hong Kong University of Science and Technology
have studied 
electrorheological fluids extensively in the context of experimentally 
exploring parameters that affect aspects of the 
transition, as well as the strength of the ER effect for 
various applications.
Figure~\ref{Ping_metal_ER}~(c)-(e) 
shows an example of their experiments, 
and the development of a fractal (DLA-like) network. 
In this experiment the electrode configuration
is a point at the bottom and a half-circle on top.
Time-lapse images from such experiments can be 
analyzed using image processing techniques
to extract particle configurations as a function of time.
This data can be analyzed to yield information
on how correlation lengths, ramification of the clusters,
connectedness, etc., depend on the field strength and other
parameters in the system. 
Building on our existing numerical tools, we can 
use the known experimental parameters to calculate the 
energy for each image and subsequently compile an energy 
landscape for the duration of the experiment. This will 
provide invaluable insight into the characterization of 
the critical energy for phase transition. Moreover, given
the evolution of these configurations, we will be able to calculate
the evolution of the spectral measure, eigenvalue spacing 
statistics, and other diagnostic data which will give insight
into our theoretical investigations into these systems.
Calculation of transport critical exponents
will be invaluable for investigating issues of universality 
of the ER fluid transition.







\begin{figure}                \centering
\includegraphics{Ping_metal_ER.eps}
   \caption{Evolution of field-induced connectedness in ER fluids. 
    Our preliminary experiments with 50 $\mu$m glass spheres are shown
      in (a) with zero field and in (b) with a 
      1000 V/mm horizontal field which induces chain formation.    
       Metal sphere experiments conducted by our
     colleagues Ping Sheng and Weijia Wen at HKUST 
      show growth of the fractal network at field
      strengths of 150 V/mm in (c), 200 V/mm in (d),
      and 250 V/mm in (e). Simulated distributions of 
        electric field strength and potential are also shown.}
      



    \label{Ping_metal_ER}
        
\end{figure}







%
%\begin{figure}               
%\centering
%\includegraphics[]{ER_fractal.eps}
%   \caption{The fractal pattern development and simulated distribution 
%             of electric field strength shown on the left and right of 
%              each panel for metal spheres in the presence of a 
%               pin-semicircle electric field with strength 
%              (a) 100 V/mm, (b) 200 V/mm, and (c) 425 V/mm.  \label{fig:ER_fractal}}
%\end{figure}

In addition to using established experimental data with 
variable particle volume fractions, table-top experiments 
will be conducted using only a few particles. 
The behavior 
of these particles in the presence of an applied electric 
field brings forth the influence of the induced dipoles and 
electrodes on each sphere in a carefully controlled 
experiment. Such a small scale study is not only simple 
and economical, but will also provide an important 
standard by which to compare numerical results. 
We will also be able to investigate in detail the metastable
states discussed above and their effect on transitional behavior.  
A preliminary experiment we conducted at the University of Utah
with dielectric glass spheres is shown in 
Figure \ref{Ping_metal_ER} (a) and (b).
Postdoctoral fellow Joyce Lin, a key participant in this
project, has extensive experimental experience working with
spheres moving in fluids. 



We also propose to extend our numerical simulations beyond 
the analysis of a given, static configuration of particles 
suspended in a medium. The time evolution of particles and 
the energy of the system under an electric field will also 
be calculated from the initial configuration. This will not 
only provide a predictive simulator for the behavior of 
electrorheological fluids, but also allow us to explore 
vast variations in the parameter space. 
Given the strong foundation of numerical techniques that 
have been developed to calculate the
local electric field in the presence of many disks or spheres, 
we will build a more comprehensive algorithm to model the 
behavior of composite materials
\cite{Klingenberg1989,Myroshnychenko2005,Cheng2008}.
Using experimental data 
from our collaborators 
as well as from experiments we will conduct, we can test and 
validate the numerics.















\section{Sea Ice Experiments and Simulations}



%characterizing spectral measures, eignevalue spacings, etc.
%collapse of the spectral gap, etc. 
%for the thermal evolution of the brine microstructure.
%There is also the interesting issue that from the tomographic
%data we can output random graphs (containing the connectivity info) 
%with tags attached to the edges indicating pore and neck sizes.
%How do we exploit this info? What does it tell us about the
%transport props. Is there a new random matrix here?


%1.) Introduce characterization of spectral measures, eigenvalue spacings etc. to measure transport properties of sea ice






We will characterize connectedness transitions 
in the microstructure 
of sea ice through experimental, analytical, and numerical 
investigations, and explore the implications for 
the critical behavior of fluid and electrical 
transport in sea ice. 
In particular, we will study the thermal 
evolution of the spectral measure and eigenvalue spacings 
of the brine microstructure. This complex microstructure, 
provided by X-ray tomography directly taken on ice cores, 
will be used to  numerically calculate the important 
spectral properties through techniques developed from 
random matrix theory, as discussed above.

While there has been progress 
\cite{Day:JPCM-96,Day:JPCM-99,Stroud:1986:PRB}
in the past on 
obtaining the spectral measure $\mu$, such as the seminal results in 
\cite{Day:JPCM-96} for the two dimensional random
resistor network,
the approach there focused on first computing the
effective conductivity $\sigma^*$.   
Then the density $\rho(z)$ of $\mu$, with 
$\mu(dz)=\rho(z)dz$, is extracted via the limit
$\rho(z)= \frac{1}{\pi}\lim_{\delta \rightarrow 0}
Im[m(z+i \delta)]$ for $z \in [0,1]$, 
and $m=\sigma^*/\sigma_2=1-F$.
Here we propose to construct the spectral measure $\mu$
directly by finding the eigenvalues $\lambda_i$ and 
eigenvectors $\psi_i$
of matrix discretizations of the operator 
$\chi_1 \Gamma \chi_1$, which can then be viewed
as a random matrix. The indicator function $\chi_1$
is determined by our tomographic data on the brine microstructure
\cite{Golden:2007:GRL,Pringle:JGR:2009}.
With $\psi_i \psi_i^{T}$ the projection
onto the $\lambda_i$ eigenspace,  
\vspace{-1ex}
\begin{equation}        \label{spectral_measure_expansion}
\mu(d \lambda)=\sum_{i=1}^N
{\kappa_i \; \delta_{\lambda_i}(d \lambda)}, \;\;\;\;
\kappa_i = 
\langle \mathbf{e}_k \cdot \psi_i \psi_i^{T} \mathbf{e}_k   \rangle,
\vspace{-1ex}
\end{equation}
%\vspace{-1ex}
where $\langle \mathbf{f} \cdot \mathbf{g} \rangle$ 
here denotes an appropriate inner product.







\begin{figure}                \centering
\includegraphics{Figure_2_brine_connectivity_new.eps}
   \caption{(a-c) X-ray CT volume renderings of the brine phase 
        within a lab-grown sea-ice single crystal 
       with S = 9.3 ppt.  The (non-collocated) 
      8 $\times$ 8 $\times$ 2 mm 
      sub-volumes in (a) and (b) illustrate changes in 
     brine connectivity 
    during warming, with no top to bottom connections  
    in (a) and one in (b). An image from
   a cylinder of height 8 mm and diameter 21 mm with many
    such connections is shown in (c). 
   A graph of the divergence of the 
    correlation length
   at the percolation threshold is shown for 
   brine in the vertical
   direction in sea ice in (d),
    with data from mapping the tomographic images
   onto random graphs \cite{Pringle:JGR:2009}.}
      
    \label{brine_micro_evol}
        
\end{figure}









We plan to develop efficient methods for computing the spectral measure
from our tomographic data. For example, we will investigate multiscale
techniques, where the underlying partial differential equation
for $\chi_1 \Gamma \chi_1$ is solved for model pore configurations, and
small-scale information is {\it up-scaled} and incorporated into 
more efficient, larger-scale discretizations.
As the tomographic data can also produce random graphs that contain information on pore size, connectivity, and channel width, we will also investigate possible discretizations of the $\chi_1 \Gamma \chi_1$ operator from the resulting adjacency matrix.

By investigating spectral measures for sea ice 
microstructure and their evolution with temperature, 
we open up the analysis of sea ice microstructure and transport
properties to sophisticated methods in 
statistical physics 
\cite{Thompson-88,Baker-1990,Golden:JMP-5627,Golden:PRL-3935} and the
theory of random matrices 
\cite{Mehta:2004:RM,Pastur:1996:AGMMP,Konig:2005:PS,Deift:1999:CPAM}.
For example, we observe a gap in the spectrum in  
Figure \ref{spectral_measures_ben} (b) for $p=0.10$.
As discussed above, 
this gap vanishes as $p \rightarrow p_c$
with a critical exponent $\Delta$, and this exponent is related
to critical exponents for transport via the classical scaling relations
of statistical mechanics \cite{Golden:PRL-3935}.
By constructing $\mu$ for sea ice directly from 
the microstructure, we can obtain
detailed information about the spectral measure, its gaps, and the 
eigenvalue spacing distribution, and their
evolution with temperature. We expect to be able to relate such
information to the critical properties of fluid, thermal, and 
electrical transport in sea ice, since the same spectral measure
arises in integral representations for these other transport
properties. 
It will be fascinating to
explore such questions for different ice types with
different threshold behavior \cite{Golden:2010:FTASI},
investigate the role of anisotropy, and the connection
of brine channels to DLA \cite{Meakin:1998}, as illustrated
in Figure \ref{erfluid_brine_dla}.



During the austral spring of 2011 aboard the Australian icebreaker
{\it Aurora Australis}, 
we will have the opportunity to travel to Antarctica and 
extract ice cores for structural analysis as described above. 
On this expedition, students and the postdoc can 
experience first-hand the importance of such 
characterizations of sea ice for microstructural and transport
transitions,
and the implications of this study for questions of climate change. 
We will likely also have opportunities to travel to the Arctic,
near Barrow, Alaska, for sea ice field experiments. 










%\vspace{1ex}









\eject





\setcounter{page}{1} 
\pagestyle{plain} 


\bibliographystyle{plain}
\bibliography{golden,cmg,opp08,murphy_joyce,furse,joyce}


\eject


\pagestyle{empty}



\begin{center}
{\Large\bf Budget Justification}
\end{center}
%\vspace{2ex}


The bulk of the funds requested are to support the participants 
in the project: summer salary for the PI Kenneth Golden, 
partial support for postdoctoral researcher Joyce Lin 
(with the possibility of partially supporting another
postdoc in the Math Department at Utah involved
in the project),
full-time support for one graduate student
(with the possibility to provide summer support for another
graduate student),
and year-round support for two undergraduate students. 

The funds for domestic 
and foreign travel are to support trips for investigators 
and students between 
the University of Utah and the Physics Department 
at Hong Kong University 
of Science and Technology to collaborate with our close
colleagues there, 
and to present results at 
domestic and 
international conferences. 
These funds will also help support the 
participation of Golden, Lin, and students in
field investigations of phase transitions in sea ice  
in the Arctic Ocean off Barrow Alaska, and in
the Antarctic sea ice pack during a scheduled Australian
voyage in 2011 aboard the icebreaker {\it Aurora Australis}.
We will also leverage other sources of funding to
help make possible these important experiences for students
in gaining a broad exposure to phase transitions in 
composites, and sea ice in particular.
We are requesting modest funding 
for consultant services to help fund short 
visits and lectures by scientists in closely related areas. 

  


We are requesting \$9,270 to cover the cost of a 
Sun Fire X4275 Server: 2 x Quad-Core Intel Xeon X5570
processor 2.93 GHz, 64GB memory, 2 x 300 GB  hard disks 
(with educational discount).
Calculations involving the
spectral measures, eigenvalue spacing densities, 
statistical mechanical functions,
and large networks 
will be highly demanding on computing resources,
and we are requesting a dedicated machine for these numerically intensive
components of the proposed work. 
We have also requested funds in materials and supplies for 
laptops, software, cameras, a power supply, a wave form generator,
glassware, and other items necessary to carry out and document
our small scale experiments on ER fluids. A portion of these funds
will be used to obtain tomographic data on the brine microstructure
of sea ice for use in our calculations on the
thermal evolution of its spectral properties.




\eject




\pagestyle{empty}




\begin{center}
 {\Large \bf Postdoctoral Mentoring Plan}
\end{center}

\noindent
The postdoctoral fellow will be based in the Math Department 
at Utah with Golden,
but will be involved in a broad range of highly interdisciplinary activities.  
These include:
theoretical and computational investigations into the 
behavior of ER fluids, sea ice, and DLA,
and the role of related sea ice transport processes 
in polar climate models; carrying out small scale experiments
on ER fluids and sea ice; collaborating with our colleagues in Physics
at the Hong Kong University of Science and Technology (including
travel to work with them there),
field work in the Arctic and Antarctic to obtain hands-on
experience with phase transitions and associated
critical behavior of transport processes in the polar
marine environment; delivering lectures in the Math Department and
at domestic and international conferences; helping to mentor REU 
and graduate students. He or she will meet formally with 
Golden once a week,
with many more informal interactions throughout the week.





\end{document}








\vspace{1ex}



\noindent
{\bf Outreach and Broader Impact:} 
As described in Results from Prior NSF Support and the Biographies, 
the investigators have been extensively involved in 
outreach to students and broader communities. We  
have been actively involved in addressing the 
broader impact of our work and its connection to 
climate change $-$ Golden is chair of the national Committee for Math
Awareness Month April 2009 on Mathematics and Climate
(www.mathaware.org), on behalf of the Joint Policy Board for Mathematics.  
Golden and Furse have both been very active in undergraduate 
education, and in attracting students into science, engineering, and math. 
We anticipate that this project will fuel even more success in this 
direction, particularly since it combines math, geophysics,
engineering, questions of climate change, and the prospect of 
Arctic and Antarctic field work.
Eicken and Golden plan to bring together mathematical and 
geophysical scientists at an AGU session on transport in sea ice. Furse
is leading a large statewide effort in Utah to bring more students into
engineering and science, and will use this project and its 
relation to climate change as a valuable tool in this effort, particularly
for student summer engineering camps, teacher training, etc. 








\vspace{0.5ex}

\noindent
{\bf Impact of the Proposed Work:} 













%compressed powder GRL, JGR




%Random resistor network

%stat mech approach

%random matrices

%DLA





% Introduce Stieltjes representations of effective parameters,
%	connections to exact representations of system energy,
%	and applications in bounds on effective parameters: 
Given the importance and long history of problems involving
transport through disordered inhomogeneous
media, going back at least to the 1800's, many
methods have been developed to solve them. 
In the late 1970's a powerful approach
was introduced to bound the effective 
transport properties of two phase 
composite media. By exploiting the analytic properties of effective parameters,
such as the complex permittivity $\epsilon^*$, this {\it analytic continuation method}
\cite{Bergman:PRC-377,Bergman:PRL-1285,Milton:RP-79,Milton:APL-300,Golden:CMP-473} 
yields a reduced representation of the system energy in terms of
a Stieltjes integral representation of $\epsilon^*$, which encapsulates all the 
complicated geometric information about the composite into a spectral measure 
$\mu$ \cite{Golden:PRL-3935}. 
The integral arises from a resolvent representation for the electric field 
involving a sparse 
self-adjoint random matrix which {\it separates} material properties from 
the applied field 
strength and geometrical information in $\mu$ \cite{Golden:PRL-3935}.

% Introduce proposed work: Statistical mechanics framework
Here we {\bf propose} to characterize phase transitions of 
connectedness and transport properties in binary composite media using
random matrix theory, orthogonal polynomial theory, and the statistical 
mechanics framework
which directly follows from the reduced energy integral representation 
\cite{Murphy_Stat_Mech}. 
Future results obtained by these analytical techniques will be motivated 
and tested by extensive numerical
simulations and small scale laboratory experiments. The elegant statistical 
mechanics framework parallels 
the Ising model and from this direct analogue many results 
have followed from standard techniques \cite{Murphy_Stat_Mech}, 
such as the relations between empirical parameters and response 
functions \cite{Baker-1990,Robertson}. 
Maxwell's relations also follow from standard techniques 
\cite{Murphy_Stat_Mech}, demonstrating that
changes in system parameters in response to other parameter variation 
is governed by rational 
functions of associated derivatives of a logarithmic functional. 
A standard perturbation analysis 
has shown that these logarithmic functionals constitute a detailed 
balance of free energy 
in the system \cite{Murphy_Stat_Mech}.  This framework is interlaced 
with many other logarithmic 
functionals from orthogonal polynomial theory \cite{Assche,Deift}, 
random matrix theory 
\cite{Deift,Mehta:2004:RM}, and large deviation theory \cite{Firas,Deift}, 
which yield valuable
analytic and geometric information in terms of their minimizers 
\cite{Murphy_Stat_Mech}. 
Based on these preliminary findings, which will be investigated 
in much greater detail, 
we plan to describe ground state configurations and phase transitions 
of connectedness 
and transport in ER fluids, percolating conductors, and sea ice in 
terms of minima and 
analytic properties of these logarithmic functionals and spectral 
transitions in $\mu$. 
 
There are many composite systems which exhibit macroscopic critical behavior as some
parameter in the system varies. For example electrorheological (ER) fluids, consisting 
of solid spherical inclusions in a dielectric liquid host, exhibit a liquid to solid phase 
transition as the applied electric field strength surpasses a critical value 
\cite{Halsey:S-761,Tao-EF,Wen:APL-2147,Wen:PF-1826}. Fractal net structures similar to 
diffusion limited aggregates form when metal inclusions are used and periodic columnar 
structures qualitatively similar to brine inclusions in sea ice form using dielectric inclusions, 
increasing the viscosity of the composite by several orders of magnitude within a few 
milliseconds leading to a solid phase \cite{Wen_etal}. If sensors are used to trigger the applied electric field, 
one can turn many devices such as clutches, valves, dampers \textit{etc}. into active 
mechanical elements capable of responding to environmental variations \cite{Wen_etal}. 
The reduced form of the system
energy, involving $\epsilon^*$, leads to an incredibly detailed description of  
energy partitions within the ER fluid in terms of Herglotz functions involving $\mu$ 
\cite{Murphy_Stat_Mech}. This yields an exact description of the evolution and exchange 
of energy as the ER fluid transitions to its solidified ground state. Preliminary numerical, 
analytical, and experimental evidence suggests the energy landscape is a hilly terrain 
associated with metastable states \cite{Murphy_Stat_Mech}, a phenomena which goes beyond the 
classical postulates of statistical mechanics \cite{Robertson}. This detailed Herglotz energy 
information describes the associated evolutions in the Gibbs--Helmholtz free energy landscape
\cite{Murphy_Stat_Mech} characterizing the critical field strength of the phase transition, 
a long standing problem, by the overcoming of local energy minima. 
 
% SEA ICE PHASE TRANSITIONS HERE 


% Introduce basic concepts of random matrix theory
Geometric resonances of two-component composites are eigenvalues of a real
symmetric sparse random matrix lying at the heart of two phase 
homogenization theory \cite{Golden:PRL-3935}. Generically, the probability 
distribution of ensembles of symmetric random matrices can be  
expressed solely in terms of its eigenvalues \cite{Mehta:2004:RM,Deift}.
Surprisingly, this distribution is in the exact form \cite{Mehta:2004:RM,Deift,Balian,Dyson}
of the statistical mechanics model of two dimensional positive charges on a 
line in an external confining potential competing with logarithmic repulsion! 
In total agreement with statistical mechanics \cite{Pastur}, the ground state distribution 
of eigenvalues is obtained by the minimum value of the associated ``electrostatic 
energy" whose distribution on the spectral interval is determined by the 
configurational minimum of the associated Gibbs-Helmholtz free energy. 
The loci of this eigenvalue distribution are also given by the large $N$ limit of 
zeros of a $N^{th}$ order polynomial \cite{Assche} which attain the unique minimum 
of a logarithmic potential, varied over all possible configurations in the support of $\mu$
\cite{Pastur}. The ground state distribution is system-specific although fluctuations 
about this distribution are universal \cite{Muttalib_etal_qRME}. These universal eigenvalue 
fluctuations couple to physical geometric fluctuations in the system which are determined 
by the eigenvalue locations and residues contained in $\mu$. The universal fluctuations in
eigenvalues of a symmetric random matrix have been exploited to obtain extremely detailed 
information about phase transitions in a wide variety of disordered systems \cite{Canali}. 


% Compare RMT mesoscopic/quantom conductors with percolating RRN
This random matrix framework has been extremely successful in 
describing the disorder--driven conduction/insulation phase transition exhibited
by mesoscopic and quantum conductors \cite{Muttalib_etal_qRME}. The phase transition
from the conducting regime to the insulating regime, with disorder increasing, is 
characterized by a transition in eigenvalue statistics from highly correlated 
Wigner-Dyson (WD) statistics, giving rise to the phenomenon of eigenvalue repulsion 
\cite{Canali}, toward uncorrelated Poisson statistics \cite{Muttalib_etal_qRME,Canali}.
Preliminary numerical simulations have demonstrated that
eigenvalue statistics of the sparse random matrix associated with 
percolating 2--D square binary bond random resistor networks have very similar features 
to the random matrix theory characterization of phase transitions in 
mesoscopic/quantum conductors \cite{Murphy_RMT}. Our results illustrate, strikingly, that 
macroscopic electric systems exhibit the same universal fluctuations as mesoscopic 
electric systems of size $10^{-6}$m and quantum systems nine orders of magnitude smaller.
Powerful analytical tools of random matrix theory such as Mehta's orthogonal polynomial method 
\cite{Mehta:2004:RM,Deift,Muttalib_etal_qRME,Canali}, supersymmetry in field theory 
\cite{Efetov}, and replica theories \cite{Kuhn} will be brought to bear on these model
systems giving us detailed information on phase transitions and ground states in two--component 
composite media.       



*************************************************


RMT

********************************************


********************************

OLD

We have found that critical behavior of binary composites can be studied using a broad range of mathematical techniques. Here we show that random matrix theory (RMT) also provides invaluable insight into this problem. A central object in homogenization theory of composites is the self adjoint random operator 
$\chi_1\Gamma\chi_1$, where $\Gamma=\nabla(-\Delta)^{-1}\nabla\cdot$ and $\chi_1$ is the characteristic function of type one medium \cite{Golden:CMP-467}. 
The action of this random integro-differential operator on a bond random resistor network (RRN) in $\mathbb{Z}^d$ is given by that of a sparse random matrix 
\cite{Golden:CMP-467}. The differential operators $\nabla$ and $\nabla\cdot$ are given by forward and backward difference operators, respectively, 
$(-\Delta)^{-1}$ is the inverse of the lattice Laplacian, and $\chi_1$ is a random diagonal matrix of ones and zeros 
\cite{Golden:CMP-467,Murphy_RMT}.  This fundamental fact allows us to incorporate powerful analytical and variational techniques of RMT into our studies of phase transitions in disordered binary composites.
 
Since the advent of RMT
this elegant framework has accurately described phase transitions in many disordered systems, examples include random graphs, quantum chaos, and disordered quantum and mesoscopic conductors \cite{Efetov,Canali,Muttalib_etal_qRME,Kuhn}. The random matrix description of all these systems have common universal features. Most well known is the transition in eigenvalues statistics from highly correlated Wigner-Dyson (WD) statistics toward uncorrelated Poisson statistics, as a function of disorder \cite{Muttalib_etal,Canali}. The strong eigenvalue correlations in WD statistics gives rise to the phenomenon of eigenvalue repulsion \cite{Mehta:2004:RM,Deift}. Our numerical simulations demonstrate eigenvalue statistics associated with the binary bond RRN on $\mathbb{Z}^2$ also has this behavior \cite{Murphy_RMT,Muttalib_etal,Evangelou:PRL:1992}(see figure). Our results illustrate, strikingly, that macroscopic electric systems exhibit the same universal fluctuations as mesoscopic electric systems of size $10^{-6}$m and quantum systems nine orders of magnitude smaller \cite{Muttalib_etal,Evangelou:PRL:1992}. Eigenvalue statistics of mesoscopic and quantum conductors are only weakly dependent on the dimension of the system \cite{Muttalib_etal,Canali}. Part of the proposed work is the generalization of our software to higher dimensions in order to determine the dependence of RRN eigenvalue statistics on the dimension, as well as calculation of other standard eigenvalue statistics \cite{Mehta:2004:RM}. We will also use analytical and variational techniques, presented below, to characterize phase transitions and ground-states in binary composite media using mesoscopic and quantum conductors as an intuitive guide. 








The probability distribution of symmetric random matrix ensembles (RMEs) can be expressed solely in terms of its eigenvalues $\{\lambda_i\}$
\cite{Mehta:2004:RM,Deift}. Surprisingly, this distribution is in the exact form \cite{Mehta:2004:RM,Deift,Balian,Dyson} of the statistical mechanics ensemble of two dimensional positive charges on a line in an external confining potential, $V>0$, competing with logarithmic repulsion: $P=Z_N^{-1}\exp(-\beta \mathcal{H}_N)$, where     
%
\begin{align}
  \mathcal{H}_N\{\lambda_1,\ldots,\lambda_{N}\}=\beta^{-1}\sum_{k=1}^{N}V(\lambda_k)-\sum_{k>l=1}^{N}\ln|\lambda_k-\lambda_l|.
\end{align}
%
Here $N$ is the size of the matrix, $Z_N$ is the partition function, and $\beta$ depends only on the symmetry class taking values 1 and 2 for real and Hermitian symmetric matrices respectively. It has been analytically proven, in the Hermitian case ($\beta=2$), that in the infinite volume limit, $N\to\infty$, when $V(\lambda)\geq\lambda$ universal WD statistics arise, which are independent of all system parameters \cite{Muttalib_etal_qRME,Canali}. At the onset of criticality $V$ violates this inequality becoming comparable to $\ln^2{\lambda}$ and the statistics deviate toward Poisson statistics, dependent only on a single disorder parameter 
\cite{Muttalib_etal_qRME,Canali}. Numerical simulations have shown that the real symmetric case ($\beta=1$) has this same behavior \cite{Canali}. In total agreement with statistical mechanics \cite{Pastur}, for 
$V(\lambda)\geq\lambda$ the ground state distribution of eigenvalues is given by the minimum "electrostatic" energy whose distribution is given by the unique equilibrium measure $\rho$ of the associated infinite volume free energy: 
$f[\rho]=\inf_m f[m]=\lim_{N\to\infty}N^{-2}\ln{Z_N}$, where 
%
\begin{align}
  \label{eq:GroundState_Measure_min}
  f[m]=
     \int_0^1\left( 
                   V(\lambda)-\frac{\beta}{2} U(\lambda;m)\right) dm(\lambda), 
\quad U(\lambda;m)=\int_0^1\ln|\lambda-\lambda^\prime|dm(\lambda),                  
\end{align}
%
varied over eigenvalue distributions on $[0,1]$. Demonstrating, once again, the ground-states of the composite are characterized by the minimum of a free energy, for this regime of disorder. We stress that the corresponding variational problem determines both the (compact) support $\sigma$ of the measure and the form of the measure \cite{Pastur}. Moreover, the corresponding Euler-Lagrange equations show that, in the infinite volume limit, the confining potential $V$ satisfies the singular integral equation 
$V^\prime(\lambda)=\int_{\sigma} d\rho(x)/(\lambda-x)$ for $\lambda\in\sigma$ \cite{Pastur}. These key properties will be aggressively pursued in our studies of phase transitions of binary composites. 

The ground-state distribution of eigenvalues depends on the form of the confining potential $V$ although fluctuations about this distribution are universal \cite{Muttalib_etal_qRME,Canali}. This fact has been exploited to obtain extremely detailed information about phase transitions in disordered systems \cite{Muttalib_etal_qRME,Canali} using powerful analytical methods. 
For Hermitian RMEs ($\beta=2$) Mehta's orthogonal polynomial method has lead to exact closed form calculation of all relevant statistics and all $n$-point correlation functions 
\cite{Mehta:2004:RM}. Thus giving exact description of phase transitions in mesoscopic and quantum conductors described by Hermitian random matrices \cite{Muttalib_etal_qRME,Canali}. Similar results are expected to follow in our real symmetric ($\beta=1$) case \cite{Muttalib_etal,Canali}. In the absence of a generalization of Mehta's method to real symmetric RMEs, we will rely on well developed super-symmetry \cite{Efetov,Canali} and replica field theory methods
\cite{Efetov,Canali,Evangelou:PRL:1992,Kuhn} describing phase transitions in eigenvalue statistics associated with real symmetric sparse random matrices.  There is a natural question regarding Pad\'e approximants of 
$F_N(s)= P_N^{[1]}(s)/P_N(s)$. If the numerator polynomial zeros, like the denominator zeros, also follow these universal eigenvalue statistics this would imply $\mu$ itself has universal features which would completely couple physical geometric fluctuations to these universal fluctuations. This is another question we hope to address with future numerical simulations. 

************************************************

We have found that critical behavior of binary composites can be studied using a broad range of mathematical techniques. Here we show that random matrix theory (RMT) also provides invaluable insight into this problem. A central object in homogenization theory of composites is the self adjoint random operator 
$\chi_1\Gamma\chi_1$, where $\Gamma=\nabla(-\Delta)^{-1}\nabla\cdot$ and $\chi_1$ is the characteristic function of type one medium \cite{Golden:CMP-467}. 
The action of this random integro-differential operator on a bond random resistor network (RRN) in $\mathbb{Z}^d$ is given by that of a sparse random matrix 
\cite{Golden:CMP-467}. The differential operators $\nabla$ and $\nabla\cdot$ are given by forward and backward difference operators, respectively, 
$(-\Delta)^{-1}$ is the inverse of the lattice Laplacian, and $\chi_1$ is a random diagonal matrix of ones and zeros 
\cite{Golden:CMP-467,Murphy_RMT}.  This fundamental fact allows us to incorporate powerful analytical and variational techniques of RMT into our studies of phase transitions in disordered binary composites.
 
Since the advent of RMT
this elegant framework has accurately described phase transitions in many disordered systems, examples include random graphs, quantum chaos, and disordered quantum and mesoscopic conductors \cite{Efetov,Canali,Muttalib_etal_qRME,Kuhn}. The random matrix description of all these systems have common universal features. Most well known is the transition in eigenvalues statistics from highly correlated Wigner-Dyson (WD) statistics toward uncorrelated Poisson statistics, as a function of disorder \cite{Muttalib_etal,Canali}. The strong eigenvalue correlations in WD statistics gives rise to the phenomenon of eigenvalue repulsion \cite{Mehta:2004:RM,Deift}. Our numerical simulations demonstrate eigenvalue statistics associated with the binary bond RRN on $\mathbb{Z}^2$ also has this behavior \cite{Murphy_RMT,Muttalib_etal,Evangelou:PRL:1992}(see figure). Our results illustrate, strikingly, that macroscopic electric systems exhibit the same universal fluctuations as mesoscopic electric systems of size $10^{-6}$m and quantum systems nine orders of magnitude smaller \cite{Muttalib_etal,Evangelou:PRL:1992}. Eigenvalue statistics of mesoscopic and quantum conductors are only weakly dependent on the dimension of the system \cite{Muttalib_etal,Canali}. Part of the proposed work is the generalization of our software to higher dimensions in order to determine the dependence of RRN eigenvalue statistics on the dimension, as well as calculation of other standard eigenvalue statistics \cite{Mehta:2004:RM}. We will also use analytical and variational techniques, presented below, to characterize phase transitions and ground-states in binary composite media using mesoscopic and quantum conductors as an intuitive guide. 

The probability distribution of symmetric random matrix ensembles (RMEs) can be expressed solely in terms of its eigenvalues $\{\lambda_i\}$
\cite{Mehta:2004:RM,Deift}. Surprisingly, this distribution is in the exact form \cite{Mehta:2004:RM,Deift,Balian,Dyson} of the statistical mechanics ensemble of two dimensional positive charges on a line in an external confining potential, $V>0$, competing with logarithmic repulsion: $P=Z_N^{-1}\exp(-\beta \mathcal{H}_N)$, where     
%
\begin{align}
  \mathcal{H}_N\{\lambda_1,\ldots,\lambda_{N}\}=\beta^{-1}\sum_{k=1}^{N}V(\lambda_k)-\sum_{k>l=1}^{N}\ln|\lambda_k-\lambda_l|.
\end{align}
%
Here $N$ is the size of the matrix, $Z_N$ is the partition function, and $\beta$ depends only on the symmetry class taking values 1 and 2 for real and Hermitian symmetric matrices respectively. It has been analytically proven, in the Hermitian case ($\beta=2$), that in the infinite volume limit, $N\to\infty$, when $V(\lambda)\geq\lambda$ universal WD statistics arise, which are independent of all system parameters \cite{Muttalib_etal_qRME,Canali}. At the onset of criticality $V$ violates this inequality becoming comparable to $\ln^2{\lambda}$ and the statistics deviate toward Poisson statistics, dependent only on a single disorder parameter 
\cite{Muttalib_etal_qRME,Canali}. Numerical simulations have shown that the real symmetric case ($\beta=1$) has this same behavior \cite{Canali}. In total agreement with statistical mechanics \cite{Pastur}, for 
$V(\lambda)\geq\lambda$ the ground state distribution of eigenvalues is given by the minimum "electrostatic" energy whose distribution is given by the unique equilibrium measure $\rho$ of the associated infinite volume free energy: 
$f[\rho]=\inf_m f[m]=\lim_{N\to\infty}N^{-2}\ln{Z_N}$, where 
%
\begin{align}
  \label{eq:GroundState_Measure_min}
  f[m]=
     \int_0^1V(\lambda)dm(\lambda)-\frac{\beta}{2}\int_0^1\int_0^1\ln|\lambda-\lambda^\prime|dm(\lambda)dm(\lambda^\prime).                   
\end{align}
%
Which again characterizes ground-states of ER fluids, for this regime of disorder, as minima of a free energy. This is a key property which we will pursue further in our studies of phase transitions of binary composites. 

The ground-state distribution of eigenvalues depends on the form of the confining potential $V$ although fluctuations about this distribution are universal \cite{Muttalib_etal_qRME,Canali}. This fact has been exploited to obtain extremely detailed information about phase transitions in disordered systems \cite{Muttalib_etal_qRME,Canali} using powerful analytical methods. 
For Hermitian RMEs ($\beta=2$) Mehta's orthogonal polynomial method has lead to exact closed form calculation of all relevant statistics and all $n$-point correlation functions 
\cite{Mehta:2004:RM}. Thus giving exact description of phase transitions in mesoscopic and quantum conductors described by Hermitian random matrices \cite{Muttalib_etal_qRME,Canali}. Similar results are expected to follow in our real symmetric ($\beta=1$) case \cite{Muttalib_etal,Canali}. In the absence of a generalization of Mehta's method to real symmetric RMEs, we will rely on well developed super-symmetry \cite{Efetov,Canali} and replica field theory methods
\cite{Efetov,Canali,Evangelou:PRL:1992,Kuhn} describing phase transitions in eigenvalue statistics associated with real symmetric sparse random matrices.  There is a natural question regarding Pad\'e approximants of 
$F_N(s)= P_N^{[1]}(s)/P_N(s)$. If the numerator polynomial zeros, like the denominator zeros, also follow these universal eigenvalue statistics this would imply $\mu$ itself has universal features which would completely couple physical geometric fluctuations to these universal fluctuations. This is another question we hope to address with future numerical simulations. 



STAT MECH



*****************************************

OLD JOYCE*********



In the statistical mechanics framework, an ER fluid system is modeled entirely by the underlying energetic contributions \cite{Robertson,Baker-1990,Ruelle}. As the geometry of an ER fluid evolves under the 
influence of an applied electric field, the associated spectral measure governs the transformation of energy \cite{Murphy_Stat_Mech} by changing to  preserve the energy constraint $\langle \mathbf{D}\cdot\mathbf{E}_f\rangle=0$ on
the system \cite{Golden:CMP-473}.
Consistent with experiment \cite{Wen_etal}, the minimum energy principle of statistical mechanics predicts that the system will evolve toward 
configurations of low intrinsic energy, due to an exponentially decreasing probability of high intrinsic energy configurations as field strength increases \cite{Robertson,Baker-1990,Ruelle}. Thus, the high field probability 
distribution is localized about the lowest energy states. Such a phenomenon is physically manifested by a 
phase transition from highly disordered (entropic) to ordered states \cite{Wen_etal}. When the electric field strength
surpasses a critical value $E_c$, structures form in the ER fluid that increase the viscosity by several orders of magnitude within a few 
milliseconds, resulting in a transition to a solid phase \cite{Wen_etal}. Generically, the equilibrium distribution of a statistical mechanics model is given 
by the global minimum of the associated Gibbs-Helmholtz free energy, and during phase transition, the analytic properties of this free energy change suddenly \cite{Robertson,Baker-1990,Ruelle}. With this mathematical foundation, we will characterize equilibrium configurations and phase 
transitions of ER fluids by the minima and analytic properties of the associated free energy.

Our 2--D square lattice dipole model of a dielectric ER fluid captures the qualitative features of 
low field chain-like states and high field columnar states \cite{Murphy_Stat_Mech} (see figure).
We will extend this numerical model by generalizing configurations to continuous systems of disks,
spheres, and ellipsoids. These simulations will obtain the energy of a configuration by direct calculation of the electric field,
providing a more accurate description of ER fluid evolutions.


%The minimum energy principle of statistical mechanics indicates why the Herglotz spectral representations of system energies are so important in the 
%modeling of ER fluids. 
In the Herglotz spectral representations of the system energies, properties of the free energy are determined by the analytic properties of an effective parameter that 
encompasses all complicated interactions, which is both mathematically tractable and physically 
transparent. The Herglotz representations of energy partitions give an exact description of the evolution and exchange of energy 
as the ER fluid transitions to ia solidified ground state. Preliminary numerical, analytical, and experimental evidence suggests 
the energy landscape during phase transition is a hilly terrain with metastable states \cite{Murphy_Stat_Mech}, which is a phenomenon going beyond the 
classical postulates of statistical mechanics \cite{Robertson}. This detailed Herglotz information regarding energy evolutions characterizes associated 
evolutions in the Gibbs--Helmholtz free energy landscape \cite{Murphy_Stat_Mech}. %We will use this detailed spectral information 
%to rigorously solve the longstanding problem of the existence of the critical field strength $E_c$. 
Using these 
sophisticated analytical methods, we will rigorously show the existence of  $E_c$, which will also be characterized as the overcoming of local energy minima with increasing applied field strength.

Toward this goal, we present here a statistical mechanics framework of two phase composite media.
We consider a finite box $B_n \subset \mathbb{R}^d$, containing $n$ particles 
in a field $\mathbf{E}_0$ and at temperature $T$.  The system Hamiltonian is given by \cite{Murphy_Stat_Mech}
%
\begin{align} 													\label{eq:ER_fluid_Hamiltonian}
\mathcal{H}_{\omega}(E_0,s)=\frac{1}{2}P_\omega^*E_0=\frac{1}{2N}\epsilon_l E_0^2\chi_{\omega}^*(s)
\end{align}
%
for any configuration $\omega\in\Omega$, where $N$ is the suspension number density of $B_n$,
$\epsilon_l$ is the permittivity of the liquid host, $P_\omega^*$ is the effective polarization, and 
$\chi_{\omega}^*(s)=F_{\omega}(s)-p_1/s$ is the effective permeability which has integral representation 
\eqref{gold1.eq.3} with $\omega$ dependent measure \cite{Murphy_Stat_Mech}. 
We are interested in the infinite volume limit of the free energy
$f(T,E_0,s)=\lim_{n\to \infty}\frac{1}{n\beta}\ln{Z_n}$, where $Z_n$ is the probability distribution normalization called the partition function, such that the particle density $N$ remains constant. Then the average polarization of the infinite system 
$P^*(T,E_0,s)=\frac{1}{N}\epsilon_l E_0\lim_{n\to\infty}\frac{1}{n}\langle\chi_{\omega}^*(s) \rangle$, where $\langle\cdot\rangle$ denotes average over 
$\omega\in\Omega$ with Gibbs weights $\exp(-\beta\mathcal{H}_\omega)$, can be expressed in terms of the free energy as $P^*(T,E_0,s)=-\partial f/\partial E_0$ \cite{Murphy_Stat_Mech}. 
The electric susceptibility $\chi^*(T,E_0,s)=\lim_{n\to\infty}\frac{1}{n}\langle\chi_{\omega}^*(s) \rangle= F(T,E_0,s)$ 
is an example of a response function. Other response functions are defined by the first law of statistical mechanics \cite{Murphy_Stat_Mech} 
%
\begin{align}
df=-SdT-P^*dE_0-\Psi ds, 
\end{align}
%
where $S$ and $\Psi$ denote the entropy and contrast response function respectively.

Maxwell's relations are obtained by equating mixed derivatives of free energy, 
e.g. $-\partial^2f/\partial sE_0=\partial P^*/\partial s=\partial\Psi/\partial E_0$ \cite{Robertson,Baker-1990}. 
An incredible consequence of Maxwell's relations is that the covariance of the
Hamiltonian and the contrast derivative is identically zero independent of details regarding the 
evolution to equilibrium, geometry class, and parameter values. 
This universal constraint on the system implies $ \text{Cov}(F_{\omega}(s),F_{\omega}^\prime(s))\equiv 0$
\cite{Murphy_Stat_Mech}, which dramatically simplifies the analysis, yields explicit 
analytical representations of response functions $S$ and $\Psi$, and connects $\Psi$ with 
the energy contained in $\langle\chi_1 E^2\rangle$ \cite{Murphy_Stat_Mech}.
Maxwell's relations also show variations in macroscopic parameters are governed by a 
logarithmic functional $\Phi(T,E_0,s)\equiv\ln{F(T,E_0,s)}$ \cite{Murphy_Stat_Mech},
e.g. 
 $ -\partial s/\partial T=\partial S/\partial\Psi=(\partial\Phi/\partial T)(\partial\Phi/\partial E_0)^{-1}$.
%\end{align}
%
These relations constitute a
detailed balance of energy in the system, which can be shown by perturbation analysis \cite{Murphy_Stat_Mech}.
This logarithmic functional takes on a more
suggestive form when we consider Pad\'{e} approximants of
$F_\omega(s)\approx P_m^{[1]}(s)/P_m(s)$.
Denoting the infinite volume root density of the polynomials $P_m$ and $P_m^{[1]}$ by
$\rho_m$ and $\rho_m^{[1]}$, respectively, we have \cite{Assche} 
%
\begin{align} \label{eq:Log_Potential}
  \Phi(T,E_0,s)=\int_0^1\ln(s-\lambda)\langle d\rho\left( \lambda \right)-d\rho^{[1]}\left( \lambda \right)\rangle
=U\left( s;\langle\rho\rangle \right)-U(s;\langle\rho^{[1]}\rangle ),
\end{align}
%
 in the infinite volume limit. Thus, well-developed techniques of logarithmic potential theory \cite{Staff_Totik,Deift} can be used to study the energetic
processes described by $\Phi$ in \eqref{eq:Log_Potential}. A key question is whether or not a minimum principle is associated with $\Phi$.
Our numerical simulations of a 2D discrete approximation have demonstrated that minimization of the logarithmic 
potential $U(s;\mu)$ yields chain-like ground states for dielectric inhomogeneities, and fractal net ground states for metal inhomogeneities where $s=1$.
This minimum principle is also of key interest in the characterization of ground-states in ER fluids.

The elegant representation of the system Hamiltonian \eqref{eq:ER_fluid_Hamiltonian} mathematically parallels the Ising model, providing guiding intuition for our future investigation of ER fluids. For example, the onset of critical behavior in the Ising model is caused by the collapse of a gap in partition function zeros. Likewise, non-analytic behavior of the free energy in binary composites is caused by the collapse of spectral gaps in $F_\omega(s)$ as $p \rightarrow p_c$ for high contrast materials, leading to non-analytic behavior in $F_\omega(s)$ \cite{Golden:PRL-3935}. 
Furthermore, just as the separation of parameters in the Ising model gives a polynomial representation of the partition function in the variable $\exp(-\beta H)$, the parameter separation in binary composites gives a fractional polynomial representation of the partition function in the variable $\exp(-\beta \epsilon_l E_0^2/2N)$. One of the fundamental questions which must be addressed is whether or not this representation yields a factorization of the partition function which leads to a logarithmic integral representation of the associated free energy, as in the Ising model. Recent work suggests the family of polynomials orthogonal to $\mu$, complete in the Hilbert space $L^2(\mu)$, may provide a natural basis to expand the partition function in order to obtain a logarithmic integral representation of the free energy. If such a representation exists, well-developed methods of logarithmic potential theory \cite{Staff_Totik,Deift} could also be used to characterize ground states of the ER fluid as minima of the potential.



************************************

When modeling ER fluids using techniques 
from statistical mechanics, it is very important to understand various underlying energies as the system is modeled 
entirely by these energetic contributions \cite{Robertson,Baker-1990,Ruelle}. As the geometry of the ER fluid evolves under the 
influence of the applied field, the associated spectral measure evolves in a way which conserves the energy constraint on
the system $\langle \mathbf{D}\cdot\mathbf{E}_f\rangle=0$ \cite{Golden:CMP-473}. Thus governing the transformation of energy from one form to another 
\cite{Murphy_Stat_Mech}.
Consistent with experiment \cite{Wen_etal}, the minimum energy principle of statistical mechanics predicts that the system will evolve toward 
configurations of low intrinsic energy, as configurations of high intrinsic energy have exponentially shrinking probabilities of 
occurrence as the applied field strength increases \cite{Robertson,Baker-1990,Ruelle}. This leads to a high field probability 
distribution which is highly localized about the lowest energy states. This is physically manifested by a 
phase transition from highly entropic (disordered) states toward highly organized states \cite{Wen_etal}. When the electric field strength
surpasses a critical value, $E_c$, structures form increasing the viscosity of the ER fluid by several orders of magnitude within a few 
milliseconds leading to a solid phase \cite{Wen_etal}. Generically, the equilibrium distribution of a statistical mechanics model is given 
by the global minimum of the associated Gibbs-Helmholtz free energy and the phase transition is mathematically manifested by a sudden change 
in the analytic properties of this free energy \cite{Robertson,Baker-1990,Ruelle}. We hope to characterize equilibrium configurations and phase 
transitions of ER fluids by minima and analytic properties of the associated free energy.

Our 2--D square lattice dipole model of a dielectric ER fluid captures the qualitative features of 
low field chain-like states and high field columnar states \cite{Murphy_Stat_Mech} (see figure).
A key part of the proposed work will be a considerable numerical generalization of this toy model to continuous systems of disks,
spheres, and ellipsoids. These simulations will obtain the energy of a configuration by direct calculation of the electric field
giving a vastly superior description of ER fluid evolutions.


The minimum energy principle of statistical mechanics indicates why the Herglotz spectral representations of system energies are so important in the 
modeling of ER fluids. The analytic 
and convexity properties of the free energy are determined by analytic properties of the effective parameter 
encompassing all the complicated interactions in an elegant way which is mathematically tractable and physically 
transparent. The Herglotz representations of energy partitions give an exact description of the evolution and exchange of energy 
as the ER fluid transitions to its solidified ground state. Preliminary numerical, analytical, and experimental evidence suggests 
the energy landscape is a hilly terrain with metastable states \cite{Murphy_Stat_Mech}, a phenomena which goes beyond the 
classical postulates of statistical mechanics \cite{Robertson}. This detailed Herglotz information regarding energy evolutions characterize associated 
evolutions in the Gibbs--Helmholtz free energy landscape \cite{Murphy_Stat_Mech}. We will use this detailed spectral information 
in hopes of rigorously solving the long standing problem of the existence of the critical field strength $E_c$. We plan to use these 
sophisticated analytical methods to characterize $E_c$ as the overcoming of local energy minima as the applied field strength increases. 

Toward this goal, we present here a statistical mechanics framework of two phase composite media.
We consider a finite box $B_n \subset \mathbb{R}^d$ containing $n$ particles 
in a field $\mathbf{E}_0$ and at temperature $T$.  The system Hamiltonian is given by \cite{Murphy_Stat_Mech}
%
\begin{align} 													\label{eq:ER_fluid_Hamiltonian}
\mathcal{H}_{\omega}(E_0,s)=\frac{1}{2}P_\omega^*E_0=\frac{1}{2N}\epsilon_l E_0^2\chi_{\omega}^*(s)
\end{align}
%
for any configuration $\omega\in\Omega$, where $N$ is the suspension number density of $B_n$,
$\epsilon_l$ is the permittivity of the liquid host, $P_\omega^*$ is the effective polarization, and 
$\chi_{\omega}^*(s)=F_{\omega}(s)-p_1/s$ is the effective permeability which has integral representation 
\eqref{gold1.eq.3} with $\omega$ dependent measure \cite{Murphy_Stat_Mech}. 
We are interested in the infinite volume limit of the free energy
$f(T,E_0,s)=\lim_{n\to \infty}\frac{1}{n\beta}\ln{Z_n}$, where $Z_n$ is the probability distribution normalization called the partition function, such that the particle density $N$ remains constant. Then the average polarization of the infinite system 
$P^*(T,E_0,s)=\frac{1}{N}\epsilon_l E_0\lim_{n\to\infty}\frac{1}{n}\langle\chi_{\omega}^*(s) \rangle$, where $\langle\cdot\rangle$ denotes average over 
$\omega\in\Omega$ with Gibbs weights $\exp(-\beta\mathcal{H}_\omega)$, can be expressed in terms of the free energy as $P^*(T,E_0,s)=-\partial f/\partial E_0$ \cite{Murphy_Stat_Mech}. 
The electric susceptibility $\chi^*(T,E_0,s)=\lim_{n\to\infty}\frac{1}{n}\langle\chi_{\omega}^*(s) \rangle= F(T,E_0,s)$ 
is an example of a response function. Other response functions are defined by the first law of statistical mechanics \cite{Murphy_Stat_Mech} 
%
\begin{align}
df=-SdT-P^*dE_0-\Psi ds, 
\end{align}
%
where $S$ and $\Psi$ denote the entropy and contrast response function respectively.

Maxwell's relations are obtained by equating mixed derivatives of free energy, 
e.g. $-\partial^2f/\partial sE_0=\partial P^*/\partial s=\partial\Psi/\partial E_0$ \cite{Robertson,Baker-1990}. 
An incredible consequence of Maxwell's relations is that the covariance of the
Hamiltonian and the contrast derivative thereof is identically zero independent of details regarding the 
evolution to equilibrium, geometry class (e.g. spherical suspensions), and parameter values. 
This universal constraint on the system implies $ \text{Cov}(F_{\omega}(s),F_{\omega}^\prime(s))\equiv 0$
\cite{Murphy_Stat_Mech}. This dramatically simplifies the analysis of the framework, yields explicit 
analytical representations of response the functions $S$ and $\Psi$, and has connected $\Psi$ with 
energy contained in $\langle\chi_1 E^2\rangle$ \cite{Murphy_Stat_Mech}.
Maxwell's relations also show variations in macroscopic parameters are governed by a 
logarithmic functional $\Phi(T,E_0,s)\equiv\ln{F(T,E_0,s)}$ \cite{Murphy_Stat_Mech},
e.g. 
 $ -\partial s/\partial T=\partial S/\partial\Psi=(\partial\Phi/\partial T)(\partial\Phi/\partial E_0)^{-1}$.
%\end{align}
%
Perturbation analysis demonstrates these relations constitute a
detailed balance of energy in the system \cite{Murphy_Stat_Mech}.
This logarithmic functional takes on a more
suggestive form when we consider Pad\'{e} approximants of
$F_\omega(s)\approx P_m^{[1]}(s)/P_m(s)$.
Denoting the infinite volume root density of the polynomials $P_m$ and $P_m^{[1]}$ by
$\rho_m$ and $\rho_m^{[1]}$, respectively, in the infinite volume limit $(m\rightarrow\infty)$ we have \cite{Assche} 
%
\begin{align} \label{eq:Log_Potential}
  \Phi(T,E_0,s)=\int_0^1\ln(s-\lambda)\langle d\rho\left( \lambda \right)-d\rho^{[1]}\left( \lambda \right)\rangle
=U\left( s;\langle\rho\rangle \right)-U(s;\langle\rho^{[1]}\rangle ).
\end{align}
%
Therefore well developed techniques of logarithmic potential theory \cite{Staff_Totik,Deift} can be used to study the energetic
processes described by $\Phi$ \eqref{eq:Log_Potential}. A key question is whether or not a minimum principle is associated with $\Phi$.
Our numerical simulations of a discrete approximation of the 2--D checkerboard have demonstrated that minimization of the logarithmic 
potential $U(s;\mu)$ yields chain like ground-states for dielectric inhomogeneities and fractal-net ground-states for metal inhomogeneities ($s=1$).
This minimum principle is also of key interest in the characterization of ground-states in ER fluids.

The elegant representation of the system Hamiltonian \eqref{eq:ER_fluid_Hamiltonian} mathematically parallels the Ising model. The intuition we gain from this parallel will guide our future investigation of ER fluids. Indeed, much like the onset of critical behavior in the Ising model is caused by the collapse of a gap in partition function zeros, in binary composites non-analytic behavior of the free energy is caused by the collapse of spectral gaps in $F_\omega(s)$, as $p \rightarrow p_c$, for high contrast materials, $s \rightarrow 0$ or 1, leading to non-analytic behavior in $F_\omega(s)$ \cite{Golden:PRL-3935}. 
Furthermore, just as the separation of parameters in the Ising model gives a polynomial representation of the partition function in the variable $\exp(-\beta H)$, the parameter separation in binary composites gives a fractional polynomial representation (a polynomial with real valued exponents) of the partition function in the variable $\exp(-\beta \epsilon_l E_0^2/2N)$. One of the fundamental questions which must be addressed is whether or not this representation yields a factorization of the partition function which leads to a logarithmic integral representation of the associated free energy, as in the Ising model. Recent work suggests the family of polynomials orthogonal to $\mu$, complete in the Hilbert space $L^2(\mu)$, may provide a natural basis to expand the partition function in order to obtain a logarithmic integral representation of the free energy. If such a representation exists, well developed methods of logarithmic potential theory \cite{Staff_Totik,Deift} could also be used to characterize groundstates of the ER fluid as minima of the potential.


caption
A simulated annealing model of 
   microstructural evolution in electrorheologial fluids.
   Long range order develops as the applied electric field is increased 
   (left to right) and exceeds a critical value 
   $-$ similar to sea ice as temperature increases. The system 
   tends to the ground state determined by the Gibbs factor
   associated with the system Hamiltonian. 



LEE-YANG (Ben)


In the Ising model, a beautiful analytic structure of the magnetization $M$ arises from the separation of the applied magnetic field strength $H$ from details of the microstructural geometry and the Lee-Yang property \cite{Lee:PR-411} of the partition function. For finite volumes $Z$ 
is a polynomial in the activity variable $z=\exp(-\beta H)$, yielding a factorization of $Z$ as a product of roots, where $\beta=1/kT$, $k$ is Boltzmann's constant, and $T$ is the absolute temperature \cite{Lee:PR-411,Ruelle}. The Lee-Yang Theorem states that the zeros of $Z$ lie on the unit circle which
leads to a logarithmic integral representation of the Gibbs-Helmholtz free energy $f=-\beta^{-1}\ln{Z}$ \cite{Baker:PRL-990}. As the temperature approaches a 
critical value, $T\rightarrow T^+_c$, a sudden change in the analytic properties of the free energy arises from the 
collapse of a gap in partition function zeros about $z=1$ $(H=0)$ \cite{Ruelle,Baker-1990}. For $T<T_c$, divergent power law behavior occurs when $H\rightarrow0$ describing the phenomena of spontaneous magnetization  \cite{Baker-1990}. 









*****************************************************

OLD
STAT MECH  NOV 14 9am



When an electric field is applied to a dielectric ER fluid, 
the  particles become polarized, creating interactive forces 
which change the underlying geometry. As the geometry evolves, 
the associated spectral measure evolves accordingly to 
preserve the energy constraint 
$\langle \mathbf{D}\cdot\mathbf{E}_f\rangle=0$ on the system. 
Consistent with experiment \cite{Wen_etal}, the minimum 
energy principle of statistical mechanics predicts that 
the system will evolve toward configurations of low 
intrinsic energy, due to an exponentially decreasing 
probability of high intrinsic energy configurations as 
field strength increases \cite{Robertson,Baker,Ruelle}. 
Thus, the high field probability distribution is localized 
about the lowest energy states. Such a phenomenon is 
physically manifested in ER fluids by a phase transition 
from highly disordered to ordered states \cite{Wen_etal}. 
When the electric field strength
surpasses a critical value $E_c$, structures form in 
the ER fluid which increase the viscosity by several 
orders of magnitude within a few 
milliseconds, resulting in a transition to a solid 
phase \cite{Wen_etal}. Generically, the equilibrium 
distribution of a statistical mechanics model is given 
by the global minimum of the associated Gibbs-Helmholtz 
free energy, and during a phase transition, the analytic 
properties of this free energy change suddenly 
\cite{Robertson,Baker,Ruelle}. 
We will characterize equilibrium 
configurations and phase 
transitions of ER fluids by the minima and analytic 
properties of the free energy.

A 2-d square lattice dipole model of a dielectric 
ER fluid that we have developed captures the qualitative features of 
low field, disordered states and high field, columnar 
states \cite{Murphy_Stat_Mech}. The evolution of the 
states of this model as the field is increased
is shown in Figure \ref{simulated_annealing}.
We will extend this discrete model 
to continuous systems of disks and
spheres. These simulations will 
obtain the energy of a configuration by direct 
calculation of the electric field,
providing a more accurate description of ER fluid 
evolutions.

%The minimum energy principle of statistical mechanics indicates why the Herglotz spectral representations of system energies are so important in the 
%modeling of ER fluids. 
By use of Herglotz spectral representations of the 
system energies, properties of the free energy are 
determined by the analytic properties of an effective 
parameter that encompasses all complicated interactions, 
which is both mathematically tractable and physically 
transparent. The Herglotz representations of energy 
partitions give an exact description of the evolution 
and exchange of energy 
as the ER fluid transitions to a solidified 
ground state. Preliminary numerical, analytical, 
and experimental evidence suggests 
the energy landscape during phase transition is a 
hilly terrain with metastable states 
\cite{Murphy_Stat_Mech}, which is a phenomenon 
going beyond the 
classical postulates of statistical 
mechanics \cite{Robertson}. This detailed Herglotz 
information regarding energy evolutions characterizes 
associated 
evolutions in the Gibbs--Helmholtz free energy 
landscape \cite{Murphy_Stat_Mech}. 
%We will use this detailed spectral information 
%to rigorously solve the longstanding problem of the existence of the critical field strength $E_c$. 
Using these 
sophisticated analytical methods, we will 
rigorously show the existence of  $E_c$, 
which will also be characterized as the 
overcoming of local energy minima with 
increasing applied field strength.

Toward this goal, we present here a statistical 
mechanics model of dielectric ER fluids which is 
extendible to the case of metallic spheres by 
letting $s\to0$.
We consider a finite box $B_n \subset \mathbb{R}^d$, 
containing $n$ spheres in an electric field 
$\mathbf{E}_0$ and at temperature $T$.  The 
system Hamiltonian is given by \cite{Murphy_Stat_Mech}
%
\begin{align} 													\label{eq:ER_fluid_Hamiltonian}
\mathcal{H}_{\omega}(E_0,s)=
\frac{1}{2}P_\omega^*E_0=
\frac{1}{2N}\epsilon_2 E_0^2\chi_{\omega}^*(s)
\end{align}
%
for any configuration $\omega\in\Omega$, where $N$ 
is the sphere number density of $B_n$,
$\epsilon_2$ is the permittivity of the liquid host, 
$P_\omega^*$ is the polarization, and 
$\chi_{\omega}^*(s)=F_{\omega}(s)-p_1/s$ is the 
electric susceptibility which has integral representation 
\eqref{gold1.eq.3} with $\omega$ dependent 
measure \cite{Murphy_Stat_Mech}. The canonical 
partition function $Z_n$ is given by
%
\begin{align}
	Z_n(T,E_0,s)=\sum_\omega 
\exp(-\beta \mathcal{H}_{\omega}(E_0,s))=
\exp(-\beta f_n(T,E_0,s)),
\end{align}
%
where $\beta=1/kT$ and $k$ is Boltzmann's constant. 
We are interested in the infinite volume limit 
of the free energy
$f(T,E_0,s)=
\lim_{n\to \infty}\frac{1}{n}f_n(T,E_0,s)$ 
such that the particle density $N$ remains 
constant. Then the average polarization of 
the infinite system 
$P^*(T,E_0,s)=\frac{1}{N}\epsilon_l E_0
\lim_{n\to\infty}\frac{1}{n}
\langle\chi_{\omega}^*(s) \rangle$, 
where $\langle\cdot\rangle$ denotes average over 
$\omega\in\Omega$ with Gibbs weights 
$\exp(-\beta\mathcal{H}_\omega)$, can be 
expressed in terms of the free energy as 
$P^*(T,E_0,s)=-\partial f/\partial E_0$ \cite{Murphy_Stat_Mech}. 
The average electric susceptibility 
$\chi^*(T,E_0,s)=\lim_{n\to\infty}
\frac{1}{n}\langle\chi_{\omega}^*(s) \rangle= F(T,E_0,s)$ 
is an example of a response function. 
Other response functions are defined by the 
first law of statistical mechanics \cite{Murphy_Stat_Mech} 
%
\begin{align}
df=-SdT-P^*dE_0-\Psi ds, 
\end{align}
%
where $S$ and $\Psi$ denote the entropy 
and contrast response function respectively.

Maxwell's relations are obtained by equating 
mixed derivatives of free energy \cite{Robertson,Baker}, e.g. 
%%
\begin{align}
	-\frac{\partial^2 f}{\partial s E_0}=
\frac{\partial P^*}{\partial s}=
\frac{\partial\Psi}{\partial E_0}. 
\end{align}
%%
Alternate representations of the effective 
permittivity \cite{Bergman:AP-78} and 
Maxwell's relations demonstrate 
the model is stable \cite{Murphy_Stat_Mech}. 
An incredible consequence of Maxwell's relations 
is that the covariance of the
Hamiltonian and the contrast derivative is 
identically zero independent of details 
regarding the 
evolution to equilibrium, geometry class, 
and parameter values. 
This universal constraint on the system 
implies 
%
\begin{align}
	\text{Cov}(F_{\omega}(s),F_{\omega}^\prime(s))\equiv 0
\end{align}
%
\cite{Murphy_Stat_Mech}, which dramatically 
simplifies the analysis in the $s$ variable, 
yields explicit 
analytical representations of response functions 
$S$ and $\Psi$, and connects $\Psi$ with 
the energy contained in $\chi_1 E^2$ \cite{Murphy_Stat_Mech}.
Maxwell's relations also demonstrate that 
variations in macroscopic parameters are governed by a 
logarithmic functional 
$\Phi(T,E_0,s)\equiv\ln{F(T,E_0,s)}$ \cite{Murphy_Stat_Mech},
e.g. 
 $ -\partial s/\partial T=\partial S/\partial\Psi
=(\partial\Phi/\partial T)(\partial\Phi/\partial E_0)^{-1}$.
%\end{align}
%
These relations constitute a
detailed balance of energy in the system, which 
can be shown by perturbation analysis \cite{Murphy_Stat_Mech}.
This logarithmic functional takes on a more
suggestive form when we consider Pad\'{e} approximants of
$F_\omega(s)\approx P_m^{[1]}(s)/P_m(s)$.
Denoting the infinite volume root density of the 
polynomials $P_m$ and $P_m^{[1]}$ by
$\rho_m$ and $\rho_m^{[1]}$, respectively, 
we have \cite{Assche} 
%
\begin{align} \label{eq:Log_Potential}
  \Phi(T,E_0,s)=\int_0^1\ln(s-\lambda)\langle 
d\rho\left( \lambda \right)-d\rho^{[1]}
\left( \lambda \right)\rangle
=U\left( s;\langle\rho\rangle \right)
-U(s;\langle\rho^{[1]}\rangle ),
\end{align}
%
in the infinite volume limit. 
Thus, well-developed techniques of logarithmic 
potential theory \cite{Staff_Totik,Deift} 
can be used to study the energetic
processes described by $\Phi$ in 
\eqref{eq:Log_Potential}. A key question is 
whether or not a minimum principle is 
associated with $\Phi$.
Our numerical simulations of a discrete 
approximation of the 2--d composite 
checkerboard have demonstrated that 
minimization of the logarithmic 
potential $U(s;\mu)$ yields chain-like 
ground states for dielectric inhomogeneities, 
and fractal net ground states for metal inhomogeneities.
This minimum principle is also of key 
interest in the characterization of ground-states 
in ER fluids.

The elegant representation of the system 
Hamiltonian \eqref{eq:ER_fluid_Hamiltonian} 
mathematically parallels the Ising model, 
providing guiding intuition for our future 
investigation of ER fluids. For example, the 
onset of critical behavior in the Ising model 
is caused by the collapse of a gap in partition 
function zeros as 
$T\to T_c$. Likewise, non-analytic behavior of 
the free energy in binary composites is caused by 
the collapse of spectral gaps in $F_\omega(s)$ as 
$p \rightarrow p_c$ for high contrast materials, 
leading to non-analytic behavior in 
$F_\omega(s)$ \cite{Golden:PRL-3935}. 
Furthermore, just as the separation of parameters 
in the Ising model gives a polynomial representation 
of the partition function in the variable 
$\exp(-\beta H)$, where $H$ is the magnetic 
field strength, the parameter separation in 
binary composites gives a fractional polynomial 
representation of the partition function in the 
variable $\exp(-\beta \epsilon_2 E_0^2/2N)$. 
One of the fundamental questions which must be 
addressed is whether or not this representation 
yields a factorization of the partition function 
which leads to a logarithmic integral representation 
of the associated free energy, as in the Ising model. 
Recent work suggests the family of polynomials 
orthogonal to $\mu$, complete in the Hilbert 
space $L^2(\mu)$, may provide a natural basis 
to expand the partition function in order to 
obtain a logarithmic integral representation 
of the free energy. If such a representation 
exists, well-developed methods of logarithmic 
potential theory \cite{Staff_Totik,Deift} could 
also be used to characterize ground states of 
the ER fluid as minima of the potential.




OLD nov 14 10 am


\section{Random Matrix Characterization of Phase Transitions Binary Composites}





\begin{figure}                \centering
\includegraphics{Figure_Unfolded_Spacing_Density_SF.eps}
   \caption{Spectral statistics of the 2-d bond RRN. (a) The nearest neighbor eigenvalue spacing distribution for various levels of volume fraction $p$. The conducting regime, $p=p_c=0.5$, is characterized by WD-like statistics which transition, with decreasing $p$, toward Poisson-like statistics in the insulating regime.  (b) The corresponding evolution of the spectral function $\mu$. For the insulating regime ($p<p_c$) gaps in the spectral function (SF) maintain analytic behavior of $\sigma^*$ even after gaps in the eigenvalue density collapse. The SF gaps subsequently collapse, with increasing $p$, giving rise to critical behavior at the percolation threshold 
$p_c$. As expected, the behavior is symmetric about $p=0.5$. For example, the figures for $p$ and $1-p$ are identical, up to a scaling factor, when the buildup at the spectral endpoints is discarded \cite{Murphy_RMT}.}  



      



    \label{spectral_measures_ben}
        
\end{figure}





We have found that critical behavior of binary composites can be studied using a broad range of mathematical techniques. Here we show that random matrix theory (RMT) also provides invaluable insight into this problem. A central object in homogenization theory of composites is the self adjoint random operator 
$\chi_1\Gamma\chi_1$, where $\Gamma=\nabla(-\Delta)^{-1}\nabla\cdot$ and $\chi_1$ is the characteristic function of type one medium \cite{Golden:CMP-467}. 
The action of this random integro-differential operator on a bond random resistor network (RRN) in $\mathbb{Z}^d$ is given by that of a sparse random matrix 
\cite{Golden:CMP-467}. The differential operators $\nabla$ and $\nabla\cdot$ are given by forward and backward difference operators, respectively, 
$(-\Delta)^{-1}$ is the inverse of the lattice Laplacian, and $\chi_1$ is a random diagonal matrix of ones and zeros 
\cite{Golden:CMP-467,Murphy_RMT}.  This fundamental fact allows us to incorporate powerful analytical and variational techniques of RMT into our studies of phase transitions in disordered binary composites.
 
Since the advent of RMT
this elegant framework has accurately described phase transitions in many disordered systems, examples include random graphs, quantum chaos, and disordered quantum and mesoscopic conductors \cite{Efetov,Canali,Muttalib_etal_qRME,Kuhn}. The random matrix description of all these systems have common universal features. Most well known is the transition in eigenvalues statistics from highly correlated Wigner-Dyson (WD) statistics toward uncorrelated Poisson statistics, as a function of disorder \cite{Muttalib_etal,Canali}. The strong eigenvalue correlations in WD statistics gives rise to the phenomenon of eigenvalue repulsion \cite{Mehta,Deift}. Our numerical simulations demonstrate eigenvalue statistics associated with the binary bond RRN on $\mathbb{Z}^2$ also has this behavior \cite{Murphy_RMT,Muttalib_etal,Evangelou:PRL:1992}(see figure). Our results illustrate, strikingly, that macroscopic electric systems exhibit the same universal fluctuations as mesoscopic electric systems of size $10^{-6}$m and quantum systems nine orders of magnitude smaller \cite{Muttalib_etal,Evangelou:PRL:1992}. Eigenvalue statistics of mesoscopic and quantum conductors are only weakly dependent on the dimension of the system \cite{Muttalib_etal,Canali}. Part of the proposed work is the generalization of our software to higher dimensions in order to determine the dependence of RRN eigenvalue statistics on the dimension, as well as calculation of other standard eigenvalue statistics \cite{Mehta}. We will also use analytical and variational techniques, presented below, to characterize phase transitions and ground-states in binary composite media using mesoscopic and quantum conductors as an intuitive guide. 






\begin{figure}                \centering
\includegraphics{DLA_Spectral_Function.eps}
   \caption{Evolution of the spectral function for a growing 
   DLA cluster. The plots exhibit a decreasing spectral gap as the cluster percolates to the boundaries. }



    \label{DLA_Spectral_Function}
        
\end{figure}








The probability distribution of symmetric random matrix ensembles (RMEs) can be expressed solely in terms of its eigenvalues $\{\lambda_i\}$
\cite{Mehta,Deift}. Surprisingly, this distribution is in the exact form \cite{Mehta,Deift,Balian,Dyson} of the statistical mechanics ensemble of two dimensional positive charges on a line in an external confining potential, $V>0$, competing with logarithmic repulsion: $P=Z_N^{-1}\exp(-\beta \mathcal{H}_N)$, where     
%
\begin{align}
  \mathcal{H}_N\{\lambda_1,\ldots,\lambda_{N}\}=\beta^{-1}\sum_{k=1}^{N}V(\lambda_k)-\sum_{k>l=1}^{N}\ln|\lambda_k-\lambda_l|.
\end{align}
%
Here $N$ is the size of the matrix, $Z_N$ is the partition function, and $\beta$ depends only on the symmetry class taking values 1 and 2 for real and Hermitian symmetric matrices respectively. It has been analytically proven, in the Hermitian case ($\beta=2$), that in the infinite volume limit, $N\to\infty$, when $V(\lambda)\geq\lambda$ universal WD statistics arise, which are independent of all system parameters \cite{Muttalib_etal_qRME,Canali}. At the onset of criticality $V$ violates this inequality becoming comparable to $\ln^2{\lambda}$ and the statistics deviate toward Poisson statistics, dependent only on a single disorder parameter 
\cite{Muttalib_etal_qRME,Canali}. Numerical simulations have shown that the real symmetric case ($\beta=1$) has this same behavior \cite{Canali}. In total agreement with statistical mechanics \cite{Pastur}, for 
$V(\lambda)\geq\lambda$ the ground state distribution of eigenvalues is given by the minimum "electrostatic" energy whose distribution is given by the unique equilibrium measure $\rho$ of the associated infinite volume free energy: 
$f[\rho]=\inf_m f[m]=\lim_{N\to\infty}N^{-2}\ln{Z_N}$, where 
%
\begin{align}
  \label{eq:GroundState_Measure_min}
  f[m]=
     \int_0^1\left( 
                   V(\lambda)-\frac{\beta}{2} U(\lambda;m)\right) dm(\lambda), 
\quad U(\lambda;m)=\int_0^1\ln|\lambda-\lambda^\prime|dm(\lambda),                  
\end{align}
%
varied over eigenvalue distributions on $[0,1]$. Demonstrating, once again, the ground-states of the composite are characterized by the minimum of a free energy, for this regime of disorder. We stress that the corresponding variational problem determines both the (compact) support $\sigma$ of the measure and the form of the measure \cite{Pastur}. Moreover, the corresponding Euler-Lagrange equations show that, in the infinite volume limit, the confining potential $V$ satisfies the singular integral equation 
$V^\prime(\lambda)=\int_{\sigma} d\rho(x)/(\lambda-x)$ for $\lambda\in\sigma$ \cite{Pastur}. These key properties will be aggressively pursued in our studies of phase transitions of binary composites. 









The ground-state distribution of eigenvalues depends on the form of the confining potential $V$ although fluctuations about this distribution are universal \cite{Muttalib_etal_qRME,Canali}. This fact has been exploited to obtain extremely detailed information about phase transitions in disordered systems \cite{Muttalib_etal_qRME,Canali} using powerful analytical methods. 
For Hermitian RMEs ($\beta=2$) Mehta's orthogonal polynomial method has lead to exact closed form calculation of all relevant statistics and all $n$-point correlation functions 
\cite{Mehta}. Thus giving exact description of phase transitions in mesoscopic and quantum conductors described by Hermitian random matrices \cite{Muttalib_etal_qRME,Canali}. Similar results are expected to follow in our real symmetric ($\beta=1$) case \cite{Muttalib_etal,Canali}. In the absence of a generalization of Mehta's method to real symmetric RMEs, we will rely on well developed super-symmetry \cite{Efetov,Canali} and replica field theory methods
\cite{Efetov,Canali,Evangelou:PRL:1992,Kuhn} describing phase transitions in eigenvalue statistics associated with real symmetric sparse random matrices.  There is a natural question regarding Pad\'e approximants of 
$F_N(s)= P_N^{[1]}(s)/P_N(s)$. If the numerator polynomial zeros, like the denominator zeros, also follow these universal eigenvalue statistics this would imply $\mu$ itself has universal features which would completely couple physical geometric fluctuations to these universal fluctuations. This is another question we hope to address with future numerical simulations. 




%where the second relation above is the analogue of (\ref{mag_free}).
%For composite materials, we shall be interested in questions of
%minimizing $\Phi$ subject to constraints on, say, some moments of
%the distribution $\mu$, or on the support of $\mu$, or on the 
%geometrical arrangement of the phases.




****************
{\bf The Golden Potential}
For a finite sample of a resistor network (or continuum model) with N resonances $\lambda_n$, or poles of $F$, we can construct a "partition function" for transport
%
\begin{align}
\mathcal{Z}_N(s)=\prod_{n=1}^N(s-\lambda_n)
\end{align}
%
whose zeros are distributed according to $\mu$ in the $s$--plane (or $\phi$ in the h-plane), which serves as the analogue of (reference to 10). 
Then we construct a finite volume "Gibbs free energy" as follows. Define the normalized counting measure of resonances by
$\nu_N=\frac{1}{N}\sum_{i=1}^N \delta_{\lambda_i}$. All the resonances are on the interval $[0,1]$, so by Helly's selection principle there is a subsequence of the probability measures $\nu_N$ that converges weakly to a probability measure $\nu$ on $[0,1]$  \cite{Assche}. Defining the "Gibbs free energy" by 
$\Phi_N(s)=-\frac{1}{N}\log{\mathcal{Z}_N(s)}=-\int_0^1 \log|s-\lambda | d\nu_N(\lambda)$. Then weak convergence implies for $s\in K\subset\mathbb{C}/[0,1]$, where 
$K$ is compact, $\lim_{N\to\infty}|\mathcal{Z}_N(s)|^{1/N}=\exp\left( \int_0^1\log|s-\lambda|d\nu(\lambda)\right)=\exp(-\Phi(s))$. It is a standard result of logarithmic potential theory that $\Phi(s)=\inf_m\left[ -\int_\sigma \log|s-\lambda|dm(\lambda) \right]$ varied over all probability measures $m$ supported on $\sigma=supp(\mu)$ \cite{Pastur,Assche}.


the counting measure
$\nu_N$ in the $s-$plane
%(or $\phi$ in the $h-$plane), 
which
serves as the analogue of (\ref{partition_function}).






**************

old intro stat mech

We now focus on the {\bf proposed research}, the questions
we plan to address, and the tools we will develop to carry out
our investigations. 
When an electric field is applied to a dielectric ER fluid, 
the particles become polarized creating interactive forces 
which change the underlying geometry. As the geometry 
evolves, the associated spectral measure evolves 
accordingly. %to preserve the energy constraint 
%$\langle \mathbf{D}\cdot\mathbf{E}_f\rangle=0$ on the system. 
Consistent with 
experiments \cite{Wen_etal}, minimum energy principles 
of statistical mechanics predict that the system will 
evolve toward configurations of low intrinsic energy, 
due to an exponentially decreasing probability of high 
intrinsic energy configurations as field strength 
increases \cite{Robertson,Ruelle,Baker}. 
%Thus, 
%the high field probability distribution is localized 
%about the lowest energy states. Such a phenomenon is 
%physically manifested in ER fluids by a phase 
When the applied field strength $E_0$
surpasses a critical value $E_c$, the system makes a 
transition from highly disordered to ordered 
states \cite{Wen_etal}. 
Generically, the equilibrium 
distribution of a statistical mechanics model is given 
by the global minimum of the associated Gibbs-Helmholtz 
free energy, and during phase transition, the 
analytic properties of this free energy change 
suddenly \cite{Robertson,Baker,Ruelle}. 
One of our main goals is to characterize 
equilibrium configurations and phase 
transitions in ER fluids by the minima and analytic 
properties of an appropriate free energy.

Stat mech intro 2

%The minimum energy principle of statistical mechanics indicates why the Herglotz spectral representations of system energies are so important in the 
%modeling of ER fluids. 
By use of Herglotz spectral representations of the 
system energies, properties of the free energy are 
determined by the analytic properties of an effective 
parameter that encompasses all complicated interactions, 
which is both mathematically tractable and physically 
transparent. The Herglotz representations of energy 
partitions give an exact description of the evolution 
and exchange of energy 
as the ER fluid transitions to a solidified ground state. 
Preliminary numerical, analytical, and experimental 
evidence suggests 
the energy landscape during phase transition is a 
hilly terrain with metastable states \cite{Murphy_Stat_Mech}, 
which is a phenomenon going beyond the 
classical postulates of statistical mechanics \cite{Robertson}. 
This detailed Herglotz information regarding energy 
evolutions characterizes associated 
evolutions in the Gibbs--Helmholtz free energy 
landscape \cite{Murphy_Stat_Mech}. 
%We will use this detailed spectral information 
%to rigorously solve the longstanding problem of 
%the existence of the critical field strength $E_c$. 
Using these sophisticated analytical methods, we will 
rigorously show the existence of  $E_c$, which will 
also be characterized as the overcoming of local energy
minima with increasing applied field strength.


