%% Copernicus Publications Manuscript Preparation Template for LaTeX Submissions
%% ---------------------------------
%% This template should be used for copernicus.cls
%% The class file and some style files are bundled in the Copernicus Latex Package, which can be downloaded from the different journal webpages.
%% For further assistance please contact Copernicus Publications at: production@copernicus.org
%% https://publications.copernicus.org/for_authors/manuscript_preparation.html


%% Please use the following documentclass and journal abbreviations for preprints and final revised papers.

%% 2-column papers and preprints
\documentclass[journal abbreviation, manuscript]{copernicus}



%% Journal abbreviations (please use the same for preprints and final revised papers)


% Advances in Geosciences (adgeo)
% Advances in Radio Science (ars)
% Advances in Science and Research (asr)
% Advances in Statistical Climatology, Meteorology and Oceanography (ascmo)
% Annales Geophysicae (angeo)
% Archives Animal Breeding (aab)
% Atmospheric Chemistry and Physics (acp)
% Atmospheric Measurement Techniques (amt)
% Biogeosciences (bg)
% Climate of the Past (cp)
% DEUQUA Special Publications (deuquasp)
% Drinking Water Engineering and Science (dwes)
% Earth Surface Dynamics (esurf)
% Earth System Dynamics (esd)
% Earth System Science Data (essd)
% E&G Quaternary Science Journal (egqsj)
% EGUsphere (egusphere) | This is only for EGUsphere preprints submitted without relation to an EGU journal.
% European Journal of Mineralogy (ejm)
% Fossil Record (fr)
% Geochronology (gchron)
% Geographica Helvetica (gh)
% Geoscience Communication (gc)
% Geoscientific Instrumentation, Methods and Data Systems (gi)
% Geoscientific Model Development (gmd)
% History of Geo- and Space Sciences (hgss)
% Hydrology and Earth System Sciences (hess)
% Journal of Bone and Joint Infection (jbji)
% Journal of Micropalaeontology (jm)
% Journal of Sensors and Sensor Systems (jsss)
% Magnetic Resonance (mr)
% Mechanical Sciences (ms)
% Natural Hazards and Earth System Sciences (nhess)
% Nonlinear Processes in Geophysics (npg)
% Ocean Science (os)
% Polarforschung - Journal of the German Society for Polar Research (polf)
% Primate Biology (pb)
% Proceedings of the International Association of Hydrological Sciences (piahs)
% Safety of Nuclear Waste Disposal (sand)
% Scientific Drilling (sd)
% SOIL (soil)
% Solid Earth (se)
% State of the Planet (sp)
% The Cryosphere (tc)
% Weather and Climate Dynamics (wcd)
% Web Ecology (we)
% Wind Energy Science (wes)


%% \usepackage commands included in the copernicus.cls:
%\usepackage[german, english]{babel}
%\usepackage{tabularx}
%\usepackage{cancel}
%\usepackage{multirow}
%\usepackage{supertabular}
%\usepackage{algorithmic}
%\usepackage{algorithm}
%\usepackage{amsthm}
%\usepackage{float}
%\usepackage{subfig}
%\usepackage{rotating}

%%%% nbenjaminmurphyOUR MACROS %%%%
\usepackage{natbib} % needed for \citep{ref} and \citet{ref} references


\usepackage{xcolor}
\usepackage{siunitx}




\newtheorem{theorem}{Theorem}
\newtheorem{corollary}{Corollary}



\newcommand\beps{\mbox{\boldmath${\epsilon}$}}
\newcommand\bkappa{\mbox{\boldmath${\kappa}$}}
\newcommand\bmu{\mbox{\boldmath${\mu}$}}
\newcommand\bnu{\mbox{\boldmath${\nu}$}}
\newcommand\bphi{\mbox{\boldmath${\phi}$}}
\newcommand\btimes{\mbox{\boldmath${\times}$}}
\newcommand\bxi{\mbox{\boldmath${\xi}$}}

\renewcommand{\d}{\mathrm{d}}
\renewcommand{\k}{\mathrm{k}}
\newcommand\vecx{\mathbf{x}}
\newcommand\vecv{\mathbf{v}}

\newcommand{\figref}[1]{Fig.\,\ref{#1}}
\newcommand{\numfigref}[1]{\ref{#1}}

\renewcommand\d{\mathrm{d}}
\renewcommand\L{\mathrm{L}}

\newcommand{\I}{\mathrm{i}}

\newtheorem{Thm}{Theorem}[section]

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

\def\e{\epsilon}
\def\s{\sigma^*}
\def\D{\partial}
\def\->{\to}
\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}{\text{v}}
\newcommand{\uu}{\text{u}}
\newcommand{\xx}{\textbf{x}}
\newcommand{\ff}{\textbf{f}}
\newcommand{\KK}{\textbf{K}}
\newcommand{\VV}{\textbf{V}}
\newcommand{\kk}{\textbf{k}}
\newcommand{\II}{\textbf{I}}

\def\vska{\vskip 0.05cm}
\def\vskm{\vskip -0.1cm}
\def\vskb{\vskip 0.2cm}
\def\noi{\noindent}

%%%Christians Marcos%%%%%
\newcommand{\bC}{\mbox{\boldmath${C}$}}
\newcommand{\bL}{\mbox{\boldmath${L}$}}
\newcommand{\Ls}{\mbox{\boldmath${L}_s$}}
\newcommand{\bId}{\mbox{\boldmath${I}$}}

%%%%%%%%%%%%%%%%%%%
%% Ben's Macros %%%
%%%%%%%%%%%%%%%%%%%
\newcommand{\Ab}{\mathbf{A}}
\newcommand{\Ac}{\mathcal{A}}
\newcommand{\Amc}{\mathsf{A}}
\newcommand{\Am}{\mathsf{A}}
\newcommand{\As}{\mathscr{A}}
\newcommand{\Bb}{\mathbf{B}}
\newcommand{\Bm}{\mathsf{B}}
\newcommand{\Cb}{\mathbf{C}}
\newcommand{\Cm}{\mathsf{C}}
\newcommand{\Db}{\mathbf{D}}
\newcommand{\Ib}{\mathbf{}}
\newcommand{\Dmc}{\mathsf \it D}
\newcommand{\Dm}{\mathsf{D}}
\newcommand{\Ds}{\mathscr{D}}
\newcommand{\Dt}{\tilde{\mathbf{D}}}
\newcommand{\Eb}{\mathbf{E}}
\newcommand{\Ec}{\mathcal{E}}
\newcommand{\Et}{\tilde{\mathbf{E}}}
\newcommand{\Fc}{\mathcal{F}}
\newcommand{\Fs}{\mathscr{F}}
\newcommand{\Gb}{\mathbf{G}}
\newcommand{\Gm}{\mathsf{G}}
\newcommand{\Hb}{\mathbf{H}}
\newcommand{\Hc}{\mathcal{H}}
\newcommand{\Hmc}{\mathsf{H}}
\newcommand{\Hm}{\mathsf{H}}
\newcommand{\Hs}{\mathscr{H}}
\newcommand{\Imag}{\mbox{Im}\,} % cf plain TeX's \Im
\newcommand{\Jb}{\mathbf{J}}
\newcommand{\Kbc}{\mbox{\boldmath${\mathcal{K}}$}}
\newcommand{\Kb}{\mathbf{K}}
\newcommand{\Kc}{\mathcal{K}}
\newcommand{\Km}{\mathsf{K}}
\newcommand{\Mb}{\mathbf{M}}
\newcommand{\Mc}{\mathcal{M}}
\newcommand{\Mm}{\mathsf{M}}
\newcommand{\Nb}{\mathbf{N}}
\newcommand{\Nc}{\mathcal{N}}
\newcommand{\Ob}{\mathbf{0}}
\newcommand{\Om}{\mathsf{O}}
\newcommand{\Pen}{\mbox{\textit{Pe}}}  % Peclet number
\newcommand{\Pm}{\mathsf{P}}
\newcommand{\Qb}{\mathbf{Q}}
\newcommand{\Qm}{\mathsf{Q}}
\newcommand{\Rb}{\mathbf{R}}
\newcommand{\Real}{\mbox{Re}\,} % cf plain TeX's \Re and Reynolds number
\newcommand{\Rm}{\mathsf{R}}
\newcommand{\Sb}{\mathbf{S}}
\newcommand{\Sc}{\mathcal{S}}
\newcommand{\Sigc}{\Sigma_{\text{cont}}}
\newcommand{\Sigp}{\Sigma_{\text{pp}}}
\newcommand{\Smc}{\mathsf{S}}
\newcommand{\Sm}{\mathsf{S}}
\newcommand{\Ss}{\mathscr{S}}
\newcommand{\Tb}{\mathbf{T}}
\newcommand{\Tc}{\mathcal{T}}
\newcommand{\Ub}{\mathbf{U}}
\newcommand{\Uc}{\mathcal{U}}
\newcommand{\Um}{\mathsf{U}}
\newcommand{\Vb}{\mathbf{V}}
\newcommand{\Vc}{\mathcal{V}}
\newcommand{\Vm}{\mathsf{V}}
\newcommand{\Wb}{\mathbf{W}}
\newcommand{\Wm}{\mathsf{W}}
\newcommand{\Wt}{\tilde{\mathbf{W}}}
\newcommand{\Xb}{\mathbf{X}}
\newcommand{\Zb}{\mathbf{Z}}
\newcommand{\Zm}{\mathsf{Z}}
\newcommand{\bDelta}{\mbox{\boldmath${\Delta}$}}
\newcommand{\bGamma}{\mbox{\boldmath${\Gamma}$}}
\newcommand{\bI}{\mathbf{I}}
\newcommand{\bLambda}{\mbox{\boldmath${\Lambda}$}}
\newcommand{\bUpsilon}{\mbox{\boldmath${\Upsilon}$}}
\newcommand{\balpha}{\mbox{\boldmath${\alpha}$}}
%\newcommand{\bcdot}{\mbox{\boldmath${\cdot}$}}
\newcommand{\bcolon}{\boldsymbol{:}}
%\newcommand{\btimes}{\mbox{\boldmath${\times}$}}
\newcommand{\bepsilon}{\mbox{\boldmath${\epsilon}$}}
\newcommand{\bpsi}{\mbox{\boldmath${\psi}$}}
\newcommand{\brho}{\mbox{\boldmath${\rho}$}}
\newcommand{\bsig}{\mbox{\boldmath${\sigma}$}}
%\newcommand{\bsig}{\mbox{\boldmath${\sigma}$}}
%\newcommand{\bxi}{\mbox{\boldmath${\xi}$}}
\newcommand{\chit}{\tilde{\chi}}
\newcommand{\corref}[1]{Corollary~\ref{#1}}
\newcommand{\gt}{\tilde{g}}
\newcommand{\lemref}[1]{Lemma~\ref{#1}}
\newcommand{\secref}[1]{Section~\ref{#1}}
\newcommand{\numsecref}[1]{\ref{#1}}
\newcommand{\numappref}[1]{\ref{#1}}
\newcommand{\appref}[1]{Appendix~\ref{#1}}
\newcommand{\thmref}[1]{Theorem~\ref{#1}} 
\newcommand{\veca}{\boldsymbol{a}}
\newcommand{\vecb}{\boldsymbol{b}}
\newcommand{\vecchi}{\boldsymbol{\chi}}
\newcommand{\vecg}{\boldsymbol{g}}
\newcommand{\veck}{\boldsymbol{k}}
\newcommand{\vecpsi}{\boldsymbol{\psi}}
\newcommand{\vecq}{\boldsymbol{q}}
\newcommand{\vecr}{\boldsymbol{r}}
\newcommand{\vecvarphi}{\boldsymbol{\varphi}}
\newcommand{\vecw}{\boldsymbol{w}}
\newcommand{\vecxi}{\boldsymbol{\xi}}
\newcommand{\vecy}{\boldsymbol{y}}
\newcommand{\veczeta}{\boldsymbol{\zeta}}
\newcommand{\vecz}{\boldsymbol{z}}
\renewcommand{\d}{\mathrm{d}}
\newcommand{\cddot}{\boldsymbol{:}}



%%%%%%%%%%%%%%%%%
%% vector font %%
%%%%%%%%%%%%%%%%%
%% bold
% \newcommand\vecE{\mathbf{E}}
% \newcommand\vecJ{\mathbf{J}}
% \newcommand\vecD{\mathbf{D}}
% \newcommand\vece{\mathbf{e}}
% \newcommand\hate{\hat{\mathbf{e}}}
% \newcommand\bnabla{\mbox{\boldmath${\nabla}$}}
% \newcommand\bcdot{\mbox{\boldmath${\cdot}$}}
% \newcommand{\vecu}{\boldsymbol{u}}


%% not bold
\newcommand\vecE{E}
\newcommand\vecJ{J}
\newcommand\vecD{D}
\newcommand\vece{e}
\renewcommand\vecw{\text{w}}
\newcommand\vecwt{\text{\emph{w}}}
\renewcommand\vecx{x}
\newcommand\hate{\hat{e}}
\newcommand\bnabla{\mbox{$\nabla$}}
\newcommand\bcdot{\mbox{${\cdot}$}}
\newcommand\vecu{u}


\begin{document}

\title{Stieltjes Functions and Spectral Analysis \\
	in the Physics of Sea Ice}


% \Author[affil]{given_name}{surname}

\Author[1]{Kenneth M.}{Golden}
\Author[1]{N. Benjamin}{Murphy}
\Author[1]{Daniel}{Hallman}
\Author[1]{Elena}{Cherkaev}

\affil[1]{Department of Mathematics, University of Utah \\ 
	155 S 1400 E RM 233, Salt Lake City, UT 84112-0090 USA}
%\affil[]{ADDRESS}

%\title{Stieltjes Functions and Spectral Analysis \\
%	in the Physics of Sea Ice}
%\author{Kenneth M. Golden}   
%\author{N. Benjamin Murphy}   
%\author{Elena Cherkaev}
%\affil{Department of Mathematics, University of Utah \\ 
%	155 S 1400 E RM 233, Salt Lake City, UT 84112-0090}

%% The [] brackets identify the author with the corresponding affiliation. 1, 2, 3, etc. should be inserted.

%% If an author is deceased, please mark the respective author name(s) with a dagger, e.g. "\Author[2,$\dag$]{Anton}{Smith}", and add a further "\affil[$\dag$]{deceased, 1 July 2019}".

%% If authors contributed equally, please mark the respective author names with an asterisk, e.g. "\Author[2,*]{Anton}{Smith}" and "\Author[3,*]{Bradley}{Miller}" and add a further affiliation: "\affil[*]{These authors contributed equally to this work.}".


\correspondence{Kenneth M. Golden (golden@math.utah.edu)}

\runningtitle{Stieltjes Functions and Spectral Analysis in the Physics of Sea Ice}

\runningauthor{Golden, Murphy, Hallman, Cherkaev}





\received{}
\pubdiscuss{} %% only important for two-stage journals
\revised{}
\accepted{}
\published{}

%% These dates will be inserted by Copernicus Publications during the typesetting process.


\firstpage{1}

\maketitle



\begin{abstract}
Polar sea ice is a critical component of Earth’s climate system. 
As a material, it is a multiscale composite of pure ice with temperature dependent
millimeter-scale brine inclusions, 
%whose volume fraction, connectedness, and geometry depend strongly on temperature, 
and centimeter-scale polycrystalline microstructure which is
largely determined by how the ice was formed. The surface layer of the 
polar oceans can be viewed as a granular composite 
of ice floes in a 
sea water host, with floe sizes ranging from centimeters 
to tens of kilometers. A principal challenge in 
modeling sea ice and its role in climate is 
%the forward homogenization problem of
how to 
use information on smaller scale structure to find the 
effective or homogenized properties on larger scales 
relevant to process studies and coarse-grained climate 
models. That is, how do you predict macroscopic 
behavior from microscopic laws, like in 
statistical mechanics and solid state physics? 
Also of great interest in climate science 
is the inverse problem of recovering 
parameters controlling small scale processes from large scale 
observations.  % like in remote sensing. 
Motivated
by sea ice remote sensing, % in the 1990s, 
the analytic continuation method 
%a powerful mathematical method 
for obtaining rigorous bounds on %integral representations for 
the homogenized coefficients of two-phase composites 
was applied to 
%finding rigorous bounds on the
the complex permittivity of sea ice, which
%The method treats 
%the bulk transport coefficients %of composites 
is a Stieltjes function of the ratio of the 
permittivities of ice and brine. 
%such as the complex permittivities of ice and brine. 
Integral representations for the 
effective parameters distill the complexities 
of the composite microgeometry into the 
spectral properties of a self-adjoint operator like 
the Hamiltonian in quantum physics.
These techniques have been extended 
to polycrystalline materials, advection diffusion processes,
and ocean waves in the sea ice cover. 
Here we discuss this powerful approach in 
%mathematical 
homogenization, % and tour its extensions,
highlighting the spectral representations and 
resolvent structure 
of the fields that are shared
by the two component theory and its extensions.
Spectral analysis of sea ice structures 
leads to a random matrix theory picture of 
percolation processes in composites, establishing parallels 
to Anderson localization and semiconductor physics, 
 providing new insights into the physics of sea ice.
\end{abstract}


%\copyrightstatement{TEXT} %% This section is optional and can be used for copyright transfers.


\introduction  %% \introduction[modified heading if necessary]
\begin{figure}                \centering
	\includegraphics[width=6.5in]{multiscale_composite_sea_ice_NPG_Oct_2022.eps}
	\caption{\textbf{Sea ice as a multiscale composite material.} 
		From left to right: millimeter-scale brine inclusions
		that form the porous microstructure of sea ice \citep{Golden:GRL:L16501};
		%\citep{Weeks:FLSI:69}; 
		centimeter-scale
		polycrystalline structure of sea ice \citep{Arcone:JGR-14281}; 
		melt ponds on Arctic sea ice in late spring and summer 
		(D. Perovich) turn the surface into a two-phase composite of
		ice and melt water;  %, U. Bergen); 
		the sea ice pack as a granular composite viewed from space (NASA),
		with ``grains'' ranging in horizontal extent from meters
		to tens of kilometers;
		the Arctic Ocean %with its sea ice cover 
		viewed from space (NASA).}
	%       Images courtesy of Don Perovich.  }

\label{multiscale_structure_sea_ice}

\end{figure}

% ONR Homog Sea Ice 2017




%The precipitous loss of nearly half the extent of the summer 
%Arctic sea ice cover over the past four decades or so,
%since satellite observations started in 1979,
%is perhaps one of the most visible
%large-scale changes on Earth’s surface connected to planetary 
%arming, with significant implications for the Arctic
%region 
%and beyond

The precipitous loss of nearly half the extent of the summer Arctic sea ice cover in recent decades is perhaps the most dramatic, large-scale development on Earth's surface that has been observed to be connected to planetary warming \citep{Stroeve2007,Stroeve_GRL_2012,Maslanik:GRL-2007,%
Notz_CMIP6,Notz_Science}.
While the response of the
sea ice pack surrounding the Antarctic continent 
to the changing climate 
has perhaps not been as clear as in the Arctic,
this past year the summer sea ice
extent set a record low \citep{Turner_GRL_2022},
followed by a new record low in February 2023.
The emerging dynamics of Earth's polar marine
environments are complex and highly variable, and hence challenging to understand and predict. Yet these challenges must be faced; the changing sea ice pack is a key component of the greater climate system and directly impacts expanding human activities in these regions. Sea ice has bearing on almost any study
of the physics or biology of the 
polar marine system, as well as on almost 
any maritime operations or logistics.
Advancing our ability to analyze, model, and
predict the behavior of sea ice is critical to improving 
projections of climate change and the response of polar 
ecosystems, and in meeting the challenges of increased human
activities in the Arctic
\citep{Golden_NAMS_2020}.







One of the fascinating, yet formidable aspects of modeling 
sea ice and its role in global climate is the sheer range
of relevant length scales $-$ over ten orders
of magnitude, from the sub-millimeter 
scale to thousands of kilometers, as indicated in 
\figref{multiscale_structure_sea_ice}.
%In Figures \ref{multiscale_structure_A} and \ref{multiscale_structure_B},
%we show principal examples of sea ice structure illustrating
%such a range of scales. 
Modeling the macroscopic behavior of
sea ice on scales appropriate for 
climate models and process studies 
%even down to the centimeter scale,
depends on understanding the properties of 
sea ice on finer scales, down to individual floes and even 
the scale of the brine inclusions which control
many of the distinct
physical characteristics of sea ice as a material.
Climate models challenge the most powerful 
supercomputers to their fullest capacity. However, 
even the largest computers still limit the %horizontal 
resolution to tens of kilometers and typically require clever 
approximations and parameterizations 
to incorporate the basic physics of sea ice
\citep{Golden_NAMS_2020,Golden:PCAM:2015,Golden:NAMS:2009}.
One of the
fundamental challenges in modeling sea ice—and a central 
theme in what follows—is how to account
for the influence of the microscale on macroscopic behavior, that is,
how to rigorously use information about smaller scales 
to predict effective behavior on larger scales.
Here we consider three
different homogenization problems
in the physics of sea ice:
the classic two-phase problem of brine inclusions 
in an ice host, sea ice as a polycrystalline material,
and advection diffusion processes such as thermal
conduction or nutrient diffusion
in the presence of, e.g. convective brine flow.
All of these questions
are also of particular interest in polar 
microbial ecology \citep{Thomas:2003:SI,Reimer_EL_2022}.


We observe that this central problem of finding 
the effective properties of sea ice
%and its role in climate
is analogous to the main focus
of statistical mechanics,
where the dynamics of molecular interactions are used to 
find the  macroscopic behavior of physical systems
\citep{Thompson-88,Christensen:2005}.
%\textcolor{red}{From Daniel: I'd like to rewrite this next sentence but I'm having trouble understanding what it's referring to:} 
Moreover, %it   %also shares fundamental similarities
%is essentially the same problem in several cases
%as the main 
mathematical homogenization theory for 
differential equations with random coefficients
%composites where 
similarly seeks to find the large scale effective behavior 
from some knowledge about the local coefficients 
\citep{MILTON:2002:TC,Torquato:RHM-02,Bensoussan:AAPS-78,%
Papanicolaou:RF-835,Kozlov:RMS-91}.
These fields of physics and applied mathematics  
provide a natural
framework for treating sea ice in
predictive models of climate,
and improving projections of how Earth's polar
ice packs may evolve in the future.


The analytic continuation method (ACM)
\citep{Bergman:PRL-1285,Milton:APL-300,%
Golden:CMP-473,Golden:IMA-97,MILTON:2002:TC},
in particular, 
yields powerful integral representations for the 
effective or homogenized transport coefficients 
of two component \citep{Golden:CMP-473} or multicomponent 
\citep{Golden:JSP-655,Golden:JMPS-333} media.
The method exploits the properties of these 
coefficients as analytic functions of ratios of
the constituent parameters for two-phase media,
such as the ratio of the electrical or thermal
conductivities, or the complex permittivities.
The geometry of the composite
microstructure is encoded into 
a self-adjoint operator $G$
through the characteristic function %\textcolor{red}{$\chi$} 
which
takes the values 1 in one component (brine) 
and 0 in the other (ice).
The key step in obtaining the integral representation,
say in the case of electrical conductivity, 
is to derive a formula for the local electric field
in terms of the resolvent of $G$, 
and then apply the spectral theorem in an appropriate Hilbert space.
This representation for the effective
conductivity (or effective complex permittivity)
achieves a complete separation between the 
component parameters in the variable, and the
geometry of the microstructure embedded in the
\textit{spectral measure} of $G$. 
%the principal mathematical object in the integral.
In a discrete model of a
composite, the operator $G$ 
becomes a random matrix,
whose eigenvalues and eigenvectors can be 
used to compute the spectral measure 
\citep{Murphy:2015:CMS:13:4:825}.


The Stieltjes or Herglotz structure of the
effective parameters and their
integral representations can be used 
%use knowledge of 
to find rigorous bounds on the 
homogenized transport coefficients 
\citep{Bergman:PRL-1285,Milton:APL-300,Golden:CMP-473,Golden:JMPS-333,%
Baker:1996:Book:Pade,MILTON:2002:TC}. 
The bounds are based on knowledge of the moments of the spectral measure
or the correlation functions of the composite microstructure.
Bounds on the complex permittivity 
of sea ice as a two-phase composite were first obtained
in the context of remote sensing 
and the mathematical analysis of sea ice electromagnetic properties
\citep{Golden:JGR-13699,Golden:TGRS-1655,Golden:TGRS-1675}.
For example, the mass of the spectral measure
is the brine volume fraction. If this is known, then one can
obtain \textit{elementary bounds} in the complex case, which reduce
to the classical arithmetic and harmonic mean bounds
for real parameters. If the microgeometry is further assumed to be 
statistically isotropic, then tighter Hashin-Shtrikman bounds
can be obtained. Tighter bounds can be obtained
when the composite is assumed to have \textit{matrix-particle} structure,
e.g., separated brine inclusions in a pure ice host
\citep{Bruno:PRSLA-353,Golden:IMA-97}. Such structure leads to
gaps in the spectrum of $G$ and tighter constraints 
on the support of the spectral measure.

In remote sensing applications of
the inverse homogenization problem 
\citep{Cherkaev:WRM-437,Cherkaev:IP-1203}, 
one uses measurements of
bulk electromagnetic behavior, such as
the effective complex permittivity, to 
obtain information about the spectral measure 
\citep{Cherkaev:IP-1203} and 
bounds on the microstructural characteristics,
such as the brine volume fraction, connectivity,
\citep{Cherkaev:1996,Cherkaev:WRM-437,Golden:TGRS-1675,Gully:PB:2007,Orum:PRSLA:2012,Cherkaev:JBiomech-345}, 
and crystal orientation \citep{Gully:PRSA:2015}. 
%and connectivity \citep{Orum:PRSLA:2012}.
The \textit{microscale} structure, which determines
the spectral measure and the homogenized coefficient, 
is thus linked to the \textit{macroscopic} behavior
via the operator $G$ and its spectral characteristics,
and \textit{vice versa}.
In the multicomponent case with three or more constituents,
the homogenized transport
coefficients are analytic functions of two or more 
complex variables, and a polydisc representation formula
was found to obtain rigorous bounds
\citep{Golden:JSP-655,Golden:JMPS-333}.

The first area of application where the ACM was
extended beyond the classical case of two component 
and multiphase composites is
diffusive transport in the presence
of a flow field, 
which is widely encountered throughout science and engineering
\citep{McLaughlin:SIAM_JAM:780,Biferale:PF:2725,Fannjiang:1994:SIAM_JAM:333,Fannjiang:1997:1033,Pavliotis:PHD_Thesis,Majda:Kramer:1999:book,Majda:1994:10.1088,Xin:2009:Fronts:9780387876832}.
In addition to thermal, saline, and nutrient transport through 
the porous microstructure of sea ice, large scale transport of ice floes
and heat are also advection diffusion processes.
%Avellaneda and Majda 
\cite{Avellaneda:PRL-753,Avellaneda:CMP-339}
found a Stieltjes integral representation for the 
effective diffusivity as a function of the P\'{e}clet number
for diffusion in an incompressible 
velocity field. Based on the approach in \citep{Golden:CMP-473},
they set up a Hilbert space framework 
and applied the spectral theorem to a resolvent representation 
involving analogues of $G$ and the electric field, where the
spectral measure depends on the geometry of the 
velocity field, and knowledge of its moments yields bounds 
on the effective diffusivity.
In \citep{Murphy:ADSTPF-2017,Murphy:ADSRF-2019} we
proved novel versions of the Stieltjes 
%integral representation 
formulas. We also
developed a framework to numerically compute the spectral measures
and a systematic method to find its moments $-$ and thus a hierarchy
of bounds \citep{Bergman:AP-78,Golden:JMPS-333} 
for both the time dependent and independent cases.

In another extension of the ACM to a large class of media,
a Stieltjes integral representation and rigorous
bounds for the effective complex permittivity of 
polycrystalline media were developed 
in \citep{Gully:PRSA:2015}, based
on a resolvent formula for the electric field and
earlier observations in 
\citep{Milton:JAP-5286,Bergman:SSP-147,MILTON:2002:TC}.
The bounds assume knowledge of the average crystal orientation and 
the complex permittivity tensor of an individual crystal grain.
In sea ice, finding the complex permittivity tensor of
an individual crystal involves homogenizing
the smaller scale brine microstructure \citep{Gully:PRSA:2015}.
The polycrystalline structure of sea ice, as characterized by
the statistics of grain size, shape, and orientation, is influenced by the 
conditions under which the ice was grown
\citep{Weeks:CM-82,Petrich:SI:2009,Untersteiner:1986-GSI}.
For example, while sea ice grown in quiescent
conditions tends to have rather large-grained
\textit{columnar} structure, when grown in more turbulent
or wavy conditions it typically has a fine-grained 
\textit{granular} structure. These distinctly
different ice types have quite different fluid
flow properties
\citep{Golden:S-2238}. %,Golden:PTG:2021}. 
Additionally, when there is a well-defined
current direction during formation, crystal orientations 
tend to be statistically anisotropic within
the horizontal plane \citep{Weeks:JGR-1137}.
This can significantly affect the sea ice
radar signature used in
measurements of sea ice thickness and other properties
used to validate climate models
\citep{Golden:JGR-8107}. %,McLean:poly:2019}.

The interaction of ocean surface waves 
with polar sea ice is a critical process in Earth's climate system; 
its accurate representation is of great importance 
for developing efficient climate models. 
Declining summer sea ice has increased wave activity and the importance of ice-ocean interactions in the Arctic \citep{Waseda:2018:8:4489}.
The Arctic marginal ice zone (MIZ), the transitional region
between dense pack ice and open ocean
characterized by strong wave-ice and
atmosphere-ice-ocean interactions,
has widened significantly %as the climate has warmed 
\citep{Strong2013c}.
These recent changes can have complex implications 
for both sea ice formation 
and melting \citep{jmse9040365}. 
Indeed, the propagation of surface waves through 
Earth's sea ice covers is a complex phenomenon 
that drives their growth and decay. 
One of the main approaches to studying
waves in sea ice which is 
valid when wavelengths are much greater than floe sizes,
is to %assume that 
model the surface layer of the ice-covered ocean
as a continuum with effective properties
\citep{Bates1980,Keller1998,Wang2010,Mosig2015}.
%based on homogenization theory for composite media, 
%including an effective viscous dissipation.
Recently,
%obtained homogenized partial differential equations for 
this fundamental 
problem in sea ice physics was homogenized
with a Stieltjes integral representation
for the effective complex viscoelasticity of the surface layer,
based on a resolvent formula for the local strain field
\citep{Sampson:PHD:2017}.
%of the  as an ice-water composite. 
The integral involves 
a spectral measure of a self-adjoint operator 
which depends on the geometry of the floe configurations.
The mass of the spectral measure 
is the area fraction of the ocean covered by sea ice,
which is a standard satellite data product known as 
the \textit{sea ice concentration} field.
If the measure's mass is known, rigorous bounds can be obtained on the 
complex viscoelastic parameter. 
%\citep{Sampson:BCV:2021}.
Previously, this effective parameter had only 
been fitted to wave data. 
%We will leave any detailed
%discussion of waves in sea ice to other publications.


Early on in our work in extending the ACM to the
above problems in sea ice physics, it was clear that the classical
approach based on bounding effective parameters using the moments
of the spectral measure would, in many cases, have limited 
effectiveness. Bounds with only a moment or two known
can be quite wide, particularly for a high
contrast in the properties of the constituents, like in sea ice. 
We then developed a framework in the classic
two-phase case for computing the spectral
measure through discretization of the relevant microstructures
and finding the eigenvalues and eigenvectors of the 
matrix representation of $G$.
By developing the mathematical foundation
for these computations \citep{Murphy:2015:CMS:13:4:825}
and studying the properties of computed spectral measures
for a broad range of sea ice and other microstructures,
like human bone
\citep{Golden:JBM:337},  %,Murphy:2016},
we discovered that eigenvalue statistics
displayed fascinating
behavior depending on the connectedness of one of the phases.

The statistical behavior of the spectrum is related to the extent that
the eigenfunctions overlap. A key example is
the Anderson theory of the metal-insulator 
transition (MIT) \citep{Anderson:PR:1958,Evers:RMP:2008}, which
provides a powerful theoretical framework for understanding 
%the conditions that determine
when a disordered medium allows electronic transport, 
and when it does not. Indeed, for large enough disorder
the electrons are localized in different places with
uncorrelated energy levels
described by Poisson statistics \citep{Shklovskii:PRB:1993,Kravtsov:PRL:1913}.
For small disorder, the wave functions are extended and overlap, 
giving rise to correlated Wigner-Dyson (WD)
statistics \citep{Shklovskii:PRB:1993,Kravtsov:PRL:1913} 
with strong level repulsion \citep{Guhr:1998:PR.299.4.189}.
In work on %based on %the Stieltjes representation for 
the effective complex permittivity 
for electromagnetic wave propagation through
two-phase composites in the long wavelength regime
(or other transport coefficients 
such as thermal or electrical conductivity),
%like sea ice, 
we found an Anderson transition in spectral characteristics
as the microstructure developed long range order 
in the approach to a 
percolation threshold \citep{Murphy:PRL:2017}. We observed transitions in localization characteristics of the 
field vectors and associated transitions in spectral behavior
from uncorrelated Poissonian statistics to universal 
(repulsive) Wigner-Dyson statistics, %in our case 
reminiscent of the Gaussian Orthogonal Ensemble (GOE)
in random matrix theory. Moreover, mobility edges appear in a manner akin to Anderson localization, where such edges mark the characteristic energies of the quantum MIT \citep{Guhr:1998:PR.299.4.189}. In \citep{Morison:2022:quasiperiodic} a novel class 
of two-phase media was introduced
$-$ \textit{twisted bilayer composites}
based on Moir\'{e} patterns,
that display exotic effective properties and dramatic
transitions in spectral behavior with very small changes in 
system parameters.


We have laid the groundwork
for rigorous mathematical modeling of
sea ice processes using
Stieltjes integral representations for homogenized parameters
in several contexts of importance 
%including its role in the climate system and in 
%remote sensing applications
\citep{Murphy:JMP:063506,Murphy:2015:CMS:13:4:825,%
Gully:PRSA:2015,Murphy:PRL:2017,Murphy:ADSTPF-2017,%
Murphy:ADSRF-2019,%Kraitzman_thermal_cond_2023,
Golden_NAMS_2020,Golden:PCAM:2015,%
Golden:NAMS:2009,Golden:IMA-97,Golden:TGRS-1655,Golden:TGRS-1675}.
We also mention a significant recent advance in obtaining
a Stieltjes integral representation for the fluid permeability
of a porous medium {\citep{bi_ou_zhang_2023}}, and an excellent, 
recent review of 
Stieltjes integrals in materials science \citep{Luger2022,Ou2022}.
The permeability result is relevant
for sea ice modeling \citep{Golden:S-2238,Golden:GRL:L16501,Golden:NAMS:2009}
and has eluded mathematical inquiry for quite some time.

We have focused here on the central role that the composite
``microgeometry'' plays $-$ via the operator $G$ and its 
spectral measure (and analogues) $-$ in determining effective behavior on scales
relevant to coarse-grained climate models and studies of sea ice
processes. The \textit{geometry} represents different composite structures
in different contexts. At the finest scales, the 
millimeter-scale brine inclusion microstructure determines the properties of individual sea ice crystallites, whose size and orientation statistics make up the 
centimeter-scale polycrystalline microgeometry.
%At the finest scales the composite geometry 
%is determined by the brine inclusion microstructure, which in turn determines
%the properties of individual sea ice crystallites, whose size and orientation
%statistics determine the polycrystalline microgeometry. 
Convective fluid flow fields help transport heat, salt, and 
nutrients, where the flow field geometry 
on centimeter to meter scales
plays the role of the composite microstructure in advection-diffusion processes.
Ponds on the surface
of melting Arctic sea ice floes define the microgeometry of the surface composite of
melt water and snow on meter to kilometer scales. 
The surface layer of the 
ocean is a composite of sea water and sea ice, with microgeometry defined by the concentration, geometry, and arrangement of the 
ice floes on scales of meters to tens of kilometers. 
Large scale ice pack dynamics and transport 
on the scale of the Arctic Ocean
are determined primarily by advective and thermal forcing.
The ``microstructure'' of these advective wind and current
fields, as well as the temperature field, can be on scales from
meters to hundreds of kilometers.

%The homogenized material 
%properties given by the 
%Stieltjes integral representations can be used 
%for a wide range of applications in sea ice physics
%and ecology, as well as in large scale global climate
%and process models.
Stieltjes integral representations provide information on a wide range of transport parameters, including electrical and thermal conductivity, complex permittivity, diffusion coefficient, and fluid permeability. These parameters can be used as direct inputs into physical, biogeochemical, and ecological models of sea ice processes and in large scale numerical models. %Information on transport parameters 
%provided by the integral representations
%that depend on the
%brine and polycrystalline microstructures, such as 
%electrical and thermal conductivity, diffusion coefficient,
%and fluid permeability can be used directly as inputs 
%into physical, biogeochemical, and ecological models of
%sea ice processes on centimeter and meter scales, as well as
%in large scale numerical models.
The interplay between homogenization 
techniques like the analytic continuation method here and models of 
phase transitions in statistical
physics \citep{Banwell_NRP_2023}
is particularly interesting across the full range of scales.
From the millimeter-scale brine inclusions \citep{Golden:S-2238,Golden:GRL:L16501},
to meter-scale melt ponds \citep{Ma:NJP:2019} and the thermal properties
of the ice pack itself, our Stieltjes representations provide
rigorous theories of how effective parameters depend on 
the constituent parameters and mixture geometries.
%\textcolor{red}{For example, the ACM can provide the homogenized thermal conductivity of the  upper layer of the ocean, a key parameter in the dynamics of the MIZ} 
%For example, 
%in a model on the scale of the Arctic Ocean
%for the dynamics of the marginal ice zone,
%the transitional region between dense pack ice and open ocean,
%a key parameter is the homogenized thermal conductivity of the 
%upper layer of the ocean as a two phase composite of ice floes
%and sea water
%\citep{Strong_MIZ_Model_2023}. 
Finally, we note that in applications of the
ACM to wave phenomena, such as the effective complex permittivity for 
electromagnetic waves propagating through the sea ice, then the theory holds in the
quasistatic regime where the wavelength in the medium is assumed to be much longer than
the microstructural scale. Typically, EM waves in the Megahertz and low Gigahertz
frequency ranges satisfy this condition. For the effective complex viscoelasticity 
of the upper layer of the ocean, the Stieltjes representation holds again in the 
quasistatic regime where the wavelength is larger than the typical floe size.




The analytic continuation method is a powerful approach 
in homogenization that provides a robust mathematical
framework for rigorously studying effective
properties in the sea ice system.
The body of work that is discussed here
will advance our sea ice modeling capabilities 
and how sea ice is represented in global climate models, 
which will improve projections of the fate of sea ice and 
the ecosystems it supports. 
Moreover, the functions we study here in the sea ice context 
share the same mathematical properties as 
effective parameters in many other areas of science 
and engineering.  
Indeed, our work based on Stieltjes integrals has advanced the physics of materials related to twisted bilayer graphene and quasicrystals 
\citep{Morison:2022:quasiperiodic}, biomedical engineering of human bone \citep{Golden:JBM:337},
and our understanding of the physical properties of general polycrystalline media in geosciences and other fields \citep{Gully:PRSA:2015}.

%Moreover, the functions we study here in the sea ice context 
%share the same mathematical properties as 
%effective parameters in many other areas of science 
%and engineering. 
%So our work advances  
%these other areas as well, as evidenced by, e.g.,
%\citep{Morison:2022:quasiperiodic}, which is closely related
%to both twisted bilayer graphene 
%and quasicrystals in the physics of materials, 
%by \citep{Golden:JBM:337}, which concerns biomedical
%engineering of human bone,
%and \citep{Gully:PRSA:2015} which applies to general
%polycrystalline media in the geosciences and other
%fields.



\section{Percolation models.}
Connectedness of one phase in a composite material is often
the principal feature of the mixture geometry which 
determines effective behavior.
For example, if highly conducting inclusions are sparsely
distributed, forming a disconnected phase within a poorly
conducting encompassing host, then the effective conductivity will 
be poor as well. However, if there are enough conducting inclusions
so that they form connected pathways through the medium, then
the effective conductivity will be much closer to that
of the inclusions. Percolation theory 
\citep{Broadbent:PCPS-629,Stauffer-92,Grimmett-89,Bunde:FDS-91}
focuses on connectedness in disordered and inhomogeneous
media, and has provided the theoretical framework
for describing the behavior of fluid flow through sea ice 
\citep{Golden:S-2238,Golden:GRL:L16501,Golden:NAMS:2009}.



Consider the $d-$dimensional integer lattice ${\mathbb{Z}}^d$, 
and the square or cubic network of bonds joining nearest
neighbor lattice sites. In the percolation model
\citep{Broadbent:PCPS-629,Stauffer-92,Grimmett-89,Bunde:FDS-91},
we assign to each bond a conductivity $\sigma_0>0$ with probability $p$,  
meaning it is open (black),
and with probability $1-p$ we assign $\sigma_0=0$, meaning it is closed.
Two examples of lattice configurations 
are shown in \figref{statistical_physics_models}
with  $p=1/3$ in (a) and $p=2/3$ in (b).
Groups of connected open bonds are 
called {\it open clusters}. 
In this model 
there is a critical probability $p_c$,
$0<p_c<1$, the {\it percolation threshold}, at which the
average cluster size diverges and an infinite cluster
appears. For the $d=2$ bond lattice $p_c=1/2$. 
For $p<p_c$ the infinite cluster density $P_{\infty}(p)=0$,
while for $p > p_c$, $P_{\infty}(p) > 0$ and
near the threshold, 
$P_{\infty}(p) \sim (p-p_c)^{\beta}$
as $p \to p_c^+$,
where $\beta$ is a universal critical exponent.
It depends only on dimension and not on the details of the lattice.
Let $x,y \in {\mathbb{Z}}^d$ and 
$\tau(x,y)$ be the probability that $x$ and $y$
belong to the same open cluster. Then for $p < p_c$,
$\tau(x,y) \sim e^{-|x-y|/\xi(p)}$, and the correlation
length $\xi(p) \sim (p_c-p)^{-\nu}$ diverges 
with a universal critical exponent $\nu$ 
as $p \to p_c^-$.
as shown in 
\figref{statistical_physics_models} (c).




The effective conductivity $\sigma^*(p)$ of the lattice, now viewed as 
a random resistor (or conductor) network,
defined via 
Kirchoff's laws, vanishes for $p<p_c$ like $P_{\infty}(p)$ since
there are no infinite pathways,
as shown in \figref{statistical_physics_models} (e).
For  $p>p_c$, $\sigma^*(p)>0$, and near $p_c$,
$\sigma^*(p)\sim\sigma_0(p-p_c)^t, \;p \to p_c^+,$
where $t$ is the conductivity critical exponent, with 
$1 \leq t \leq 2$ in $d=2,3$ 
\citep{Golden:PRL-2923,Golden:CMP-467,Golden:HPM-27},
and numerical values $t \approx 1.3$ in $d=2$ and 
$t \approx 2.0$
in $d=3$ \citep{Stauffer-92}.
% $\sigma^*(p)$ in $d=2$ is shown in \figref{perc_fig_perm_data} (a).
Consider a random pipe network with %effective 
fluid
permeability $k^*(p)$ exhibiting similar behavior
$k^*(p)\sim k_0 (p-p_c)^e,$ 
where $e$ is the
permeability critical exponent, with $e=t$
\citep{Chayes:CMP-133,Sahimi-95,Golden:HPM-27}. Both $t$ and $e$ are believed to
be universal -- they depend only on dimension and not the lattice. 
Continuum
models such as the Swiss cheese model, 
can exhibit nonuniversal behavior with exponents different 
from the lattice case and $e \neq t$
\citep{Halperin:PRL-2391,Feng:PRB-197,Stauffer-92,
	Sahimi-94,Kerstein:JPA-3071}.







\begin{figure}                \centering
	\includegraphics[width=\textwidth]{statistical_physics_models_new.eps}
	\caption{The two dimensional square lattice percolation model below
		its percolation threshold of $p_c=1/2$ in (a) and above it in (b).
		(c) Divergence of the correlation length as $p$ approaches $p_c$.
		The infinite cluster density of the percolation model is shown in (d), and
		the effective conductivity is shown in (e). 
		% (f) The magnetization of the Ising model. 
	}
	\label{statistical_physics_models}
\end{figure}




%\vspace{2ex}
%******************************* Begin Ben's Writing
%*******************************


\section{Analytic continuation for two-phase composites.}
\label{sec:two_component}
%





We now describe the \emph{analytic continuation method} (ACM)
for studying the effective properties of composites
\citep{Bergman:PRL-1285,Milton:APL-300,Golden:CMP-473,Golden:IMA-97}.
This method has been used to obtain rigorous bounds on bulk
transport coefficients of composite materials from partial knowledge
of the microstructure, such as the volume fractions
of the phases. Examples of transport coefficients
to which this approach applies include
the complex permittivity, electrical and thermal conductivity,
diffusivity, magnetic permeability, and elasticity  \citep{Cherkaev:IP-1203,Cherkaev:JBiomech-345,Cherkaev:AAM:2019,Ou:MMAS:655,Ou:two-param,Kantor:1982:JMPS,Avellaneda:PRL-753,Avellaneda:CMP-339,Murphy:JMP:2020,Murphy:ADSTPF-2017}.
In \citep{Golden:JGR-13699,Golden:TGRS-1655,Golden:TGRS-1675,Golden:IMA-97,%
	Golden:PCAM:2015,Golden:NAMS:2009,Golden_NAMS_2020}, rigorous
bounds on the complex permittivity of sea ice were found.
%We have extended the method to 
%polycrystalline composites
%\citep{Gully:PRSA:2015,Barabash:JPCM:10323},
%\citep{Gully:PRSA:2015,Barabash:JPCM:10323} 
%advection-enhanced %\emph{bulk}
%diffusion processes
%by brine flow velocity fields
%\citep{Avellaneda:PRL-753,Avellaneda:CMP-339,Murphy:JMP:2020,Murphy:ADSTPF-2017},
%our recent result on 
%thermal transport in sea ice with
%brine convection \citep{Kraitzman_thermal_cond_2023},
%and ocean surface wave propagation through the ice cover
%\citep{Sampson:BCV:2021}.

%discussed below in Section \ref{sec:Eff_Diff} 
%as a part of our principal scientific investigations.




To set ideas we focus on the complex permittivity, 
keeping in mind the broad applicability of the ACM.
Consider a two-phase random medium with
%$\epsilon(x,\omega)$ the local permittivity,
%nbmurphy
local permittivity tensor $\beps(x,\omega)$,
a spatially stationary random field
in $x\in\mathbb{R}^d$ and $\omega\in\O$, where $\O$ is the
set of realizations of the medium. 
%Later, we will consider a polycrystalline medium where 
% the local permittivity 
%$\beps$ is a non-trivial symmetric matrix. 
We consider a two-phase locally isotropic medium, where the components $\epsilon_{jk}$,
$j,k=1,..,d$, of $\beps$
satisfy 
% $\epsilon_{jk}(x,\omega)=\epsilon(x,\omega)\,\delta_{jk}$,
%
\begin{align}\label{eq:locally_iso_sigma}
	\epsilon_{jk}(x,\omega)=\epsilon(x,\omega)\,\delta_{jk}\,,
\end{align}
%
where $d$ is dimension, $\delta_{jk}$ is the Kronecker delta, and 
%
\begin{align}\label{eq:sigma_two-phase}
	\epsilon(x,\omega) = \epsilon_1 \,\chi_1(x,\omega) 
	+ \epsilon_2 \,\chi_2(x,\omega)\,.
\end{align}
%
Later, we will consider a polycrystalline medium where $\beps$ is a non-trivial symmetric matrix. 
% Here we consider a two-phase locally isotropic medium, where 
% $\beps(x,\omega)=\epsilon(x,\omega)I$, $I$ is the identity matrix, and
% %the components 
% %$\epsilon_{jk}$ of $\beps$ satisfy $\epsilon_{jk}(x,\omega)=\epsilon(x,\omega)\,\delta_{jk}$,
% %where $\delta_{jk}$ is the Kronecker delta and 
% %with probability measure $P$.
% %$\epsilon(x) = \epsilon_1 \chi_1 + \epsilon_2 \chi_2$,
% %nbmurphy
% %$\epsilon(x,\omega) = \epsilon_1 \chi_1 + \epsilon_2 \chi_2$,
% %
% %\vspace{-1ex}
% %\begin{align}\label{eq:sigma_two-phase}
%  $ \epsilon(x,\omega) = \epsilon_1 \,\chi_1(x,\omega) 
%                   + \epsilon_2 \,\chi_2(x,\omega)$.
% %\end{align}
% %\vspace{-3.5ex}
%
Here, $\chi_i(x,\omega)$ is the characteristic function of medium $i=1,2$,
equaling 1 for $\omega \in \O$
with medium $i$ at $x$, and 0 otherwise, with $\chi_1+\chi_2=1$. 
%Let $\vecE(x)$ and $\vecD(x)$ be the stationary
%nbmurphy
The random electric and displacement fields $\vecE(x,\omega)$ and $\vecD(x,\omega)$ satisfy
% Let $\vecE(x,\omega)$ and $\vecD(x,\omega)$ be the stationary
% random electric and displacement fields, 
% satisfying
%
\begin{align}\label{eq:Maxwells_Equations}
	\bnabla\times\vecE=0\,, 
	\quad
	\bnabla\cdot\vecD = 0\,, 
	\quad
	\vecD=\beps\vecE\,.
	%   \qquad
	%   \langle\vecE\rangle = \vece_k\,,
\end{align}
%
% related by $\vecD= \epsilon \vecE$, satisfying
% %\begin{equation}                                 \label{gold1_1}  
% $\bnabla \cdot \vecD = 0$ and $\bnabla \times \vecE = 0$,                               
%\end{equation}
%where the average $\langle \vecE (x) \rangle = \vece_k$,   
A variational problem \citep{Golden:CMP-473} establishes 
that $\vecE$ can be written as $\vecE=\vecE_f+\vecE_0$ 
satisfying
%
\begin{align}\label{eq:E_variational}
	\vecE=\vecE_f+\vecE_0\,,
	\qquad
	\bnabla\times\vecE_f=0\,,    
	\quad
	\langle\vecD\cdot\vecE_f\rangle=0\,,
	\quad
	\langle\vecE\rangle=E_0\,.
\end{align}
%
In simpler terms, this says that curl-free and divergence-free fields are orthogonal (Helmholtz's theorem).
%and is rigorously established via the Lax-Milgram theorem .


The effective permittivity tensor 
$\beps^*$ is defined as
$\langle \vecD \rangle=  \beps^* \langle \vecE \rangle$,
% %
% \begin{align}\label{eq:sigma*}
	% \langle \vecD \rangle =  \beps^* \langle \vecE \rangle\,,
	% \quad
	% \langle \vecE \rangle=\vece_k\,,
	% % \langle \vecE \rangle=E_0\vece_k\,,
	% \end{align}
% %
where $\langle\cdot\rangle$ is ensemble averaging over $\Omega$ 
or, 
by an ergodic theorem,
spatial average over all of $\mathbb{R}^d$ \citep{Golden:CMP-473}.
We prescribe that $\vecE_0$ has direction $\vece_k$, the 
$k$th direction unit vector, and focus on the diagonal coefficient
$\epsilon^* = \epsilon^*_{kk}$, with
%$\epsilon^* =  \langle \epsilon E_k \rangle$, 
$\epsilon^* =  \langle \epsilon \vecE\cdot \vece_k \rangle$.
% Equivalently, it can be defined  
% in terms of system energy using $\langle \vecD\bcdot \vecE\rangle=\epsilon^*E_0^2$ \citep{Golden:CMP-473}. 
% Because of \eqref{eq:E_variational}, $\epsilon^*$ is also given in terms of the system energy \citep{Golden:CMP-473}
% %
% %\vspace{-1ex}
% \begin{align}\label{eq:sigma*}
	%   \langle \vecD\cdot \vecE \rangle =  \epsilon^*E_0^2\,,
	%   \quad
	%   \vecE_0=\langle \vecE \rangle\,,   
	%   \quad
	%   \vecE_0=E_0\vece_k\,.
	% \end{align}
% \vspace{-4ex}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Herglotz function description %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Since
% $\epsilon^*$ depends on 
% $h = \epsilon_1/\epsilon_2$, we define
% $m(h) =\epsilon^*/\epsilon_2$, which is 
% a Herglotz function. 
% It is analytic
% off $(- \infty, 0]$ in the complex $h$--plane, and 
% maps the upper half plane to the 
% upper half plane \citep{Bergman:PRC-377,Milton:APL-300,Golden:CMP-473}. 
%
%
% \noindent
% The linear constitutive
% relation $\vecD=\beps \vecE$ in \eqref{eq:Maxwells_Equations} relates the
% local intensity and flux,  while the linear relation in
% \eqref{eq:sigma*} relates the mean intensity and mean flux.
%
%
The key step 
of the method 
is to obtain the following Stieltjes integral representation for $\epsilon^*$ \citep{Bergman:PRC-377,Milton:APL-300,Golden:CMP-473,MILTON:2002:TC},  
% Consider $F(s) = 1 - m(h)$, $s = 1/(1-h)$,
% which is analytic off $[0,1]$ in the $s$--plane. Then \citep{Golden:CMP-473}
%
\begin{equation} \label{eq:F(s)}
	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 Stieltjes measure with support in $[0,1]$, and $F$ plays the role of a (negative) electric 
susceptibility \citep{Bergman:PRC-377}. 
In the variable $s=1/(1-h)$, with $h=\epsilon_1/\epsilon_2$, $F(s)$ is a \emph{Stieltjes function} \citep{Golden:PRL-3935,Cherkaev:IP-1203,Murphy:JMP:063506}.
This representation arises 
% from writing \eqref{eq:sigma_two-phase} as $\epsilon=\epsilon_2(1-1/s)$ and
%
% and the energy in \eqref{eq:sigma*} as $\langle \vecD\cdot \vecE\rangle=\epsilon_2(\langle \vecE\cdot \vecE_0\rangle+\langle \chi_1 \vecE\cdot \vecE_0\rangle/s)$ %
from
a resolvent
formula for the electric field 
%obtained via manipulation of (\ref{gold1_1}) and (\ref{gold1_2}),
%nbmurphy
(in medium 1) \citep{Murphy:2015:CMS:13:4:825,Cherkaev:IP-1203,Golden:CMP-473},
%$\chi_1 \vecE = (s - \chi_1\Gamma \chi_1)^{-1}\chi_1e_k$,
%\vspace{-0.5ex}
%
\begin{align}\label{eq:Electric_Field_two_component}
	\chi_1\vecE=s(sI-G)^{-1}\chi_1\vece_k\,,
	%\chi_1\vecE=s(sI-\chi_1\Gamma\chi_1)^{-1}\chi_1e_k\,,
	%
	% \chi_1\vecE=(sI-G)^{-1}\xi\,,
	\quad
	G=\chi_1\Gamma\chi_1,
	% \quad
	% \xi = s\chi_1e_k,
\end{align}
%
%\vspace{-0.5ex}
yielding 
%\vspace{-1ex}
%$F(s) = \langle \chi_1 [(s + \Gamma \chi_1)^{-1} \vece_k] \cdot \vece_k \rangle.$
%nbmurphy
$F(s) = \langle [(sI - G)^{-1} \chi_1 \vece_k] \cdot \vece_k \rangle$,
%\vspace{-1ex}
where $\Gamma  = -\bnabla (- \Delta )^{-1} \bnabla \cdot$ is a projection onto the range of the gradient operator $\bnabla$
and $\vece_k$ is the standard basis vector in the $k$th direction. 
% and is based on convolution with the Green's function for the Laplacian $\Delta=\nabla^2$. 
%The operator $\Gamma$ acts as a
%projection onto curl-free fields, and is based on
%convolution with the 
%free-space Green's function associated with $-\Delta$. 
% In the Hilbert space 
% $L^2 (\Omega)$, 
% %$L^2 (\Omega, P)$, 
% %with weight $\chi_1$ in the inner product, $\Gamma \chi_1$
% %nbmurphy
% the operator $G=\chi_1\Gamma \chi_1$ 
% %is bounded and 
% is self adjoint. 
% with $\parallel \Gamma \chi_1 \parallel \leq 1$.
Equation (\ref{eq:F(s)})
is the spectral representation of the resolvent formula in \eqref{eq:Electric_Field_two_component}
and $\mu$ is a spectral measure of 
%$\Gamma \chi_1$, 
the self-adjoint operator
$G=\chi_1\Gamma \chi_1$
on $L^2(\Omega,P)$.
% in the $(\vece_k,\vece_k)$ state.


A critical feature of 
equation \eqref{eq:F(s)} is that
the component parameters 
in $s$ are separated from the geometrical information in $\mu$.
% (Extensions to multicomponent media involve several complex variables
% \citep{Golden:JSP-655,Golden:JMPS-333,Milton:CPAM-647,Milton:CMP-281,DellAntonio:JMP-2688}.)
Information about the geometry enters
through the moments 
%$
%\mu_n = 
%\int^1_0 z^n d \mu (z)  =
%(-1)^n \langle \chi_1 [( \Gamma \chi_1 )^n \vece_k ] \cdot \vece_k \rangle.
%$
%
\begin{align}\label{eq:Moments_Fs}
	\mu_n = 
	\int^1_0 \lambda^n d \mu (\lambda)  =
	%\langle [\chi_1\Gamma \chi_1]^{\,n} \vece_k \cdot \vece_k \rangle,
	\langle G^{\,n} \chi_1 \vece_k \cdot \chi_1\vece_k \rangle.
	%\quad
	%G=\chi_1 \Gamma \chi_1. 
\end{align}
%
Then $\mu_0 = \phi$, where $\phi$ is the volume or area fraction
of phase 1, such as the brine volume fraction,
the open water area fraction or 
melt pond coverage 
% in Fig. \ref{spectral_measure} (d)-(f), 
and
$\mu_1 = \phi(1-\phi)/d$ if the material is statistically isotropic.  
In general, $\mu_n$ depends on the $(n+1)$--point correlation function of the
medium.
% begin add nbmurphy
This 
integral representation yields rigorous \emph{forward bounds} for the 
effective parameters of composites, given partial information on the 
microgeometry via the $\mu_n$ 
\citep{Bergman:PRL-1285,Milton:APL-300,Golden:CMP-473,Bergman:AP-78}.
The integral representations can also yield \emph{inverse bounds}, allowing one to use data about the electromagnetic response
of a sample to bound its structural parameters such as
the volume fraction of each of the components 
\citep{McPhedran:APA:19,McPhedran:MRSP:1990:195,Cherkaev:1996,Cherkaev:WRM-437,Cherkaev:IP-1203,%Cherkaev:IP-065008,
	Zhang:JCP-5390,Bonifasi-Lista:PMB-3063,Cherkaev:JBiomech-345,Day:1999,Golden:JBM:337,Cherkaev:Royal:2020}. 
See Section \ref{inverse} for more details.
% end add nbmurphy

\subsection{Spectral measure computations for two-phase composites} 
\label{sec:two_component_discrete}
%
Computing the spectral measure $\mu$ for a given 2D composite microstructure 
%geometry 
first involves discretizing a two-phase image of the composite 
into a square lattice filled with 1's and 0's corresponding to the two 
phases. On this square lattice, the action of the differential operators 
$\bnabla$ and $\bnabla\cdot$ are defined in terms of forward and backward 
difference operators \citep{Golden:CMP-467, Murphy:2015:CMS:13:4:825}. 
Then the key operator 
$\chi_1\Gamma\chi_1$, which depends on the geometry of the 
network via $\chi_1$, becomes a 
%self adjoint matrix. 
%nbmurphy
real-symmetric matrix $M$ \citep{Murphy:2015:CMS:13:4:825}.
% The spectral measure may be calculated directly from the
% eigenvalues  $\{\lambda_i\}$ and eigenvectors $\{v_i\}$ of this matrix 
% via %the formula
% % Then $\mu$ may be calculated directly from the
% % eigenvalues  $\{\lambda_i\}$ and eigenvectors $\{v_i\}$ of this matrix 
% % via \citep{Murphy:2015:CMS:13:4:825}
% % %
% % \begin{equation}                    \label{Spectral_Measure}
	% %   d\mu(\lambda)=\sum_{i=1}^{n} m_i \delta(\lambda-\lambda_i) d\lambda,\quad
	% %   m_i=\langle  e_0^T v_i v_i^T e_0  \rangle, 
	% % \end{equation}
% % %
% %nbmurphy
% %
% \begin{equation}                    \label{Spectral_Measure}
	%   d\mu(\lambda)=\sum_i \langle m_i \, 
	%                         \delta(\lambda-\lambda_i)\rangle \,
	%                         d\lambda,
	%   \quad m_i=e_0^T v_i v_i^T e_0 , 
	% \end{equation}
% %
% where $\delta(\lambda)$ is the Dirac delta function and $e_0$ is a vector of
% ones and zeros. In \citep{Golden:JB:2010} we computed 
% the spectral measure for samples of healthy and osteoporotic bone
% to distinguish them via the connectivity of the trabecular architecture.
%We also developed regularized inversion schemes to recover $\mu$
%from data on the effective complex permittivity of bone samples.
%
%
% As the 
% %system size $N$
% %nbmurphy
% network size $L$ 
% increases, the eigenvalues become
% increasingly dense in the spectral interval [0,1]. 
% The presence or absence of gaps in the spectrum 
% near the endpoints of [0,1],
% and the details of how large a gap is or how large the spectral
% values $m_i$ are, give important information pertaining to
% the connectivity and effective transport properties of the system. 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% In this discrete setting of a square lattice, 
% as in Fig. \ref{spectral_measure} (a)--(b),
% the action of the operator 
% $G$ 
% %$G=\chi_1\Gamma\chi_1$
% %$\chi_1\Gamma\chi_1$
% %the random operator $\chi_1\Gamma\chi_1$ 
% is given by that of a real-symmetric \emph{random  matrix} $M$, 
%
Here $\Gamma$ is a (non-random) projection matrix which depends only on the lattice topology and boundary conditions, and $\chi_1$ is a diagonal (random) projection matrix which determines the geometry and component connectivity of the composite medium \citep{Gully:PRSA:2015,Murphy:2015:CMS:13:4:825}. Another 
spectral approach
to finding effective properties based on analytic continuation
relies on computation of 
the electromagnetic eigenstates of individual inclusions 
\citep{Bergman:2020:PRA,Bergman:2022:PRA}.

% In \figref{fig:Gamma_Matrix_2D_3D} we visually display the matrix entries of $\Gamma$ with periodic boundary conditions. 
% for dimension $d=2$ and $d=3$. 
% The clear block matrix structure determines the properties of the components $\mu_{jk}$ of the matrix valued spectral measure $\bmu$ associated with anisotripic random media \citep{Murphy:2015:CMS:13:4:825}.



 The powerful integral representation in
 Equation (\ref{eq:F(s)}) is formulated in a continuum setting.
 However, in order to actually compute the spectral measures, we
 must discretize, for example, an image of the brine inclusions
 or ice floes onto a lattice, and represent 
 the operator $G$ as a matrix.
The following theorem provides a rigorous mathematical formulation of
integral representations for the effective parameters for finite
lattice approximations of two component media, and tells how
to compute the spectral measures. The electric
field decomposition in this theorem is established using the fundamental theorem of 
linear algebra and orthogonality properties of the ranges and kernels of matrix 
representations for $\nabla$, $\nabla\times$, and $\nabla\cdot$ \citep{Huang:2019:AMSA.4.1}
and will be published elsewhere. The integral representation
in equation \eqref{eq:Stieltjes_F_Discrete_two-component} is established in Theorem 2.1 of
\citep{Murphy:2015:CMS:13:4:825}. 
%
\begin{theorem} \label{thm:Discrete_Spectral_Theorem_ACM_two-comp}
	For each $\omega\in\Omega$, let $M(\omega)=W(\omega)\Lambda(\omega)\,W(\omega)$ be the eigenvalue
	decomposition of the real-symmetric matrix
	$M(\omega)=\chi_1(\omega)\,\Gamma\,\chi_1(\omega)$. Here, the columns of the matrix $W(\omega)$
	consist of the orthonormal eigenvectors $\vecwt_i(\omega)$, $i=1,\ldots,N$,
	of $M(\omega)$ and the diagonal matrix $\Lambda(\omega)={\rm diag}(\lambda_1(\omega),\ldots,\lambda_N(\omega))$
	involves its eigenvalues $\lambda_i(\omega)$. 
	Denote $Q_i=\vecwt_i\,\vecwt_i^{\;T}$ the projection matrix onto the 
	eigenspace spanned by $\vecwt_i$ and denote $\delta_{\lambda_i}(\d\lambda)$
	the Dirac $\delta$-measure centered at $\lambda_i$.
	The electric field $\vecE(\omega)$ satisfies $\vecE(\omega)=\vecE_0+\vecE_f(\omega)$, with
	$\vecE_0=\langle\vecE(\omega)\rangle$, $\Gamma\vecE(\omega)=\vecE_f(\omega)$, 
        and $\chi_1\vecE$ satisfies the resolvent formula in equation \eqref{eq:Electric_Field_two_component}. The
	effective complex permittivity tensor $\beps^*$ has components
	$\epsilon^*_{jk}$,
	%   $\epsilon_{jk}^*=\epsilon_2\,m_{jk}(h)$, 
	$j,k=1,\ldots,d$,  which satisfy       
	%
	\begin{align}\label{eq:Stieltjes_F_Discrete_two-component}
		&\epsilon^*_{jk}=\epsilon_2(\delta_{jk}-F_{jk}(s)), 
		&&F_{jk}(s)=\int_0^1\frac{\d\mu_{jk}(\lambda)}{s-\lambda}\,, 
		&&\d\mu_{jk}(\lambda)=\sum_{i=1}^N\langle \delta_{\lambda_i}(\d\lambda)\;\chi_1\,Q_i\hat{e}_j\cdot\hat{e}_k\rangle\,.  
	\end{align}
	% 
\end{theorem}







%  From Theorem \ref{thm:Discrete_Spectral_Theorem_ACM_two-comp}, the integral and
%  $\chi_1\vecE$ in equations \eqref{eq:F(s)} and \eqref{eq:Electric_Field_two_component} 
%  have explicit representations in terms of the eigenvalues $\lambda_i$ and eigenvectors
%  $\vecw_i$ of $M$ \citep{Murphy:2015:CMS:13:4:825},
% %
% \vspace{-1ex}
% \begin{align}\label{eq:Discrete_E_F(s)}
	% %
	% \chi_1 \vecE=s\sum_i\frac{\sqrt[]{m_i}}{s-\lambda_i}\;\vecw_i\,,
	% %\chi_1 \vecE=sE_0\sum_i\frac{\sqrt[]{m_i}}{s-\lambda_i}\;\vecw_i,
	% %
	% \qquad
	% %
	% F(s)=\sum_i\left\langle
	% \frac{m_i}{s-\lambda_i}\right\rangle\,,
	% %
	% \qquad
	%   m_i = |\vecw_i\,\cdot\,\chi_1 \hat{\vece}_k|^2,
	% %
	% \end{align}
% \vspace{-2ex}

% \noindent
% where $\hat{\vece}_k$ plays the role of a standard basis vector on the
% lattice. 
% Then the discrete spectral measure $\mu$ is given 
% %explicitly 
% in terms of the eigenvalues $\lambda_i$ and 
% orthonormal eigenvectors $\vecw_i$ of the matrix $M$
% \citep{Murphy:2015:CMS:13:4:825},   
% %  
% \vspace{-1ex}
% \begin{align}\label{eq:Stieltjes_F_Discrete_two-component}
	%   \mu(d\lambda)=\langle Q(d\lambda)\hat{\vece}_k\cdot\hat{\vece}_k\rangle,
	%   \quad 
	%   Q(d\lambda)=\sum_i\delta_{\lambda_i}(d\lambda)\chi_1Q_i\,,
	%   \quad
	%   Q_i=\vecw_i\vecw_i^T\,.
	% \end{align}
% \vspace{-2.5ex}
%
% Here,  $Q(d\lambda)$ is the 
% %(unique) 
% projection valued measure associated
% with $M$, $\delta_{\lambda_i}(d\lambda)$ is the delta measure centered at
% $\lambda_i$, and the matrix $Q_i=\vecw_i\vecw_i^T$ is a projection onto the
% eigenspace spanned by $\vecw_i$ \citep{Murphy:2015:CMS:13:4:825}.
%In the continuum setting, the 
%spectral measure of $G$ still
%has the form
%$\mu(d\lambda)=\langle Q(d\lambda)\vece_k\cdot \vece_k\rangle$ 
%there is a one-to-one correspondence between $G$ and a projection 
%valued measure $Q(d\lambda)$ 
%but $Q(d\lambda)$ 
%associated with $G$ 
%is abstract in nature and is
%not given explicitly as in \eqref{eq:Stieltjes_F_Discrete_two-component} \citep{Stone:64}. 




% \begin{figure}%[h!] 
	% \centering
	% \includegraphics[width=6.5in]{SF_and_RMT_for_RRN_and_melt_ponds.eps}
	% \caption{\textbf{Spectral measures and eigenvalue spacing 
			% statistics for the RRN and sea ice melt ponds.} 
		% As the probability $p$ of an open bond (black) in the 2D percolation model increases, configurations of open bonds become increasingly connected in (a)--(c); as the melt pond area fraction $\phi$ increases the ice/water composite becomes increasingly connected in (d)--(f) (images courtesy of Don Perovich), with corresponding spectral functions $\mu(\lambda)$ 
		% (histogram representations of the spectral measure) and eigenvalue spacing distributions (ESD) $P(z)$ below. 
		% The area under the graph of each spectral function in (d)--(f) is the melt water fraction $\phi=1-\psi$. 
		% In the bottom row we see the transition to universal Wigner-Dyson statistics as the conducting (a)--(c) and ocean (d)--(f) phases  become connected over the scale of the sample, and a delta function singularity in the spectral function $\mu(\lambda)$ develops at $\lambda=0$.}
	%   %The eigenvalues in the
	%   %ESD have been rescaled by an unfolding
	%   %procedure to have unit mean spacing. 
	%   %For the network, the width of
	%   %the gaps in the spectrum near $\lambda=0,1$ for $p=0.1$ shrink to 0 as
	%   %the percolation threshold $p=p_c=0.5$ is approached, with a
	%   %subsequent build up of spectra at $\lambda=0$ as $p$ surpasses $p_c$. The
	%   %spectral functions for the ice pack exhibit a similar behavior.}
% \label{spectral_measure}      
% \end{figure}

%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering
\includegraphics[scale=0.9]{brine_SM_E.png}
\caption{\emph{Electric fields and spectral measures for sea ice brine microstructure.
%
(a) Electric field values in log10 scale for X-ray CT images of 2D vertical cross sections of sea ice brine microstructure, with increasing connectivity from left to right. (b) Corresponding spectral functions $\mu(\lambda)$ (histogram representations) of the spectral measure $\mu$ in (c) and (d), which display spectral weights $m_i^2$ versus the associated eigenvalues $\lambda_i$ of the matrix $G$ and index $i$, respectively. The vertical bars in (d) delineate the $\delta$-functions at the spectral endpoints $\lambda=0,1$ seen in (c) from the rest of the spectrum, where $\lambda_i \lesssim 10^{-14}$ and $1-\lambda_i \lesssim 10^{-14}$. As the percentage of brine increases, the fluid phase becomes increasingly connected, resulting in a substantial increase in the strength of the electric field, with "hot spots" forming in geometric bottlenecks. Macroscopic connectivity of the brine phase is characterized by the mass of the $\delta$-function at $\lambda=0$ switching from numerically-zero, with the $m_i^2\lesssim 10^{-30}$, to $m_i^2\sim1$, giving rise to the "hot spots" in $E$ via equation \eqref{eq:Discrete_E_F(s)}. The electrical permittivity is taken to be $\varepsilon_1=63.3+i1930$ for brine and $\varepsilon_2=3.06$ for ice \citep{Backstrom:2006:CRST}. $E_0$ is taken to be vertically oriented.}}
%This provides an electrical signature of brine connectivity, with a 
%substantial increase in the strength of the electric field as the 
%system attains global connectivity. Here, $E_0$ is taken to be vertically 
%oriented.}
\label{fig:brine_SM_E}
\end{figure}



\begin{figure}
\centering
\includegraphics[width=\textwidth]{melt_ponds_thermal_fields.png}
\caption{\emph{Temperature gradient fields and spectral functions for sea ice melt pond microstructure.}
%
(a) Temperature gradient field values in log10 scale for melt pond microstructure atop Arctic sea ice, with increasing connectivity from left to right (images courtesy of Don Perovich), with corresponding spectral functions $\mu(\lambda)$ displayed below in (b). Here the gradient of temperature $T$, $-\nabla T$, plays the role of $E_f$ in equations \eqref{eq:Maxwells_Equations} and \eqref{eq:E_variational} \citep{MILTON:2002:TC} ($E_f=-\nabla\varphi$ for some electrical potential $\varphi$), the component thermal conductivities $\kappa_i$, $i=1,2$, play the role of the complex electrical permittivities $\epsilon_i$, and the heat current $Q$ plays the role of the displacement field $D$. Analogous to \figref{fig:brine_SM_E}, as the fluid phase becomes increasingly connected on macroscopic length scales, a buildup of spectral measure mass at $\lambda=0$ shown in (b) leads to the formation of a $\delta$-function at $\lambda=0$, with corresponding switching in the values of the $m_i$ from numerically-zero, with the $m_i\lesssim 10^{-30}$ for the left an middle figure panels, to $m_i\sim1$ for the rightmost panel. The $\delta$-function at $\lambda=1$ is also analogous to that in \figref{fig:brine_SM_E}. Like $E_0$ in \figref{fig:brine_SM_E}, we take the average thermal gradient to be vertically oriented. The thermal conductivity is taken to be $\kappa_1=0.5606$ \si{\watt\per\meter\per\kelvin} for melt ponds and $\kappa_2=0.3073$ \si{\watt\per\meter\per\kelvin} for the surrounding snow \citep{yen:1981:review}.  }
\label{fig:melt_ponds_SM_E}
\end{figure}



\begin{figure}
\centering
\includegraphics[width=\textwidth]{pack_ice_SM_E.png}
\caption{\emph{Temperature gradient fields and spectral functions for Arctic pack ice microstructure.}
%
(a) Temperature gradient field values in log10 scale for Arctic pack ice microstructure, with increasing connectivity from left to right (images courtesy of Don Perovich), with corresponding spectral functions $\mu(\lambda)$ displayed below in (b). Analogous to \figref{fig:brine_SM_E}, as the fluid phase becomes increasingly connected on macroscopic length scales, a buildup of spectral measure mass at $\lambda=0$ shown in (b) leads to the formation of a $\delta$-function at $\lambda=0$. We see a corresponding switching in the values of the $m_i$ from numerically-zero, with the $m_i\lesssim 10^{-30}$ for the left and middle figure panels, to $m_i\sim1$ for the rightmost panel. The $\delta$-function at $\lambda=1$ is also analogous to that in \figref{fig:brine_SM_E}. Like $E_0$ in \figref{fig:brine_SM_E}, we take the average thermal gradient to be vertically oriented.  Thermal conductivity is taken to be $\kappa_1=0.57$ \si{\watt\per\meter\per\kelvin} for ocean and $\kappa_2=2.11$ \si{\watt\per\meter\per\kelvin} for ice floes \citep{Pringle:2006:JGR:Oceans}. }
\label{fig:pack_ice_SM_E}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%   



% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% description of projection method %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
From Theorem \ref{thm:Discrete_Spectral_Theorem_ACM_two-comp}, the integral and
$\chi_1\vecE$ in equations \eqref{eq:F(s)} and \eqref{eq:Electric_Field_two_component} 
have explicit representations in terms of the eigenvalues $\lambda_i$ and eigenvectors
$\vecw_i$ of $M$ \citep{Murphy:2015:CMS:13:4:825,Murphy:PRL:2017},
%
% \vspace{-1ex}
\begin{align}\label{eq:Discrete_E_F(s)}
%
\chi_1 \vecE=s\sum_i\frac{m_i}{s-\lambda_i}\;\vecw_i\,,
%\chi_1 \vecE=sE_0\sum_i\frac{\sqrt[]{m_i}}{s-\lambda_i}\;\vecw_i,
%
\qquad
%
F(s)=\sum_i\left\langle
\frac{m_i^2}{s-\lambda_i}\right\rangle\,,
%
\qquad
m_i = \chi_1\vecw_i\,\cdot\hat{\vece}_k,
%
\end{align}
where $\hat{\vece}_k$ plays the role of a standard basis vector on the
lattice. 
To 
%numerically 
compute $\mu$, a non-standard generalization of the spectral theorem for matrices is required, due to the projective nature of the matrices $\chi_1$ and $\Gamma$ \citep{Murphy:2015:CMS:13:4:825}. 
% In particular, 
We developed a \emph{projection method} that
% demonstrates 
shows
the 
% %discrete 
spectral measure $\mu$ in \eqref{eq:Stieltjes_F_Discrete_two-component} depends only on the eigenvalues and eigenvectors of random sub-matrices of $\Gamma$ of size $N_1\approx \phi N$. They correspond to diagonal components $[\chi_1]_{ii}=1$, as the spectral weights $m_i$ (Christoffel numbers) associated with eigenvectors satisfying $\chi_1\vecw_i=0$ are themselves zero, $m_i=0$ \citep{Murphy:2015:CMS:13:4:825}. Fortunately, since these submatrices are much smaller for low volume fractions, this method greatly improves the efficiency and accuracy of numerical computations of $\mu$.
%the eigenvalues and eigenvectors of a smaller submatrix of the whole
%, which greatly improves the efficiency and accuracy of numerical %computations of $\mu$  



The measure $\mu$ exhibits fascinating transitional behavior as a
function of system connectivity.
% , which controls critical behavior of $\sigma^*$ near connectedness thresholds.
For example, in the case
of a random resistor network (RRN)
%random resistor network (RRN) 
with a low volume fraction $p$ of open bonds, as shown in
\figref{statistical_physics_models}a, there are spectrum-free
regions at the spectral endpoints $\lambda=0,1$
\citep{Murphy:JMP:063506,Murphy:2015:CMS:13:4:825}. However, as $p$ approaches the 
percolation threshold $p_c$~\citep{Stauffer-92,Torquato:RHM-02}
and the system becomes increasingly connected, these spectral gaps shrink and
then vanish~\citep{Murphy:JMP:063506,Jonckheere_Luck_JPA_1998}, 
% as shown in \figref{spectral_measure}a--c, 
leading
to the formation of $\delta$-components of $\mu$ at the spectral endpoints,
\emph{precisely}~\citep{Murphy:JMP:063506} when $p=p_c$ (and
$p=1-p_c$ in $d=3$). This leads to critical behavior of $\sigma^*$ for
insulating/conducting and conducting/superconducting
systems \citep{Murphy:JMP:063506}. 
This gap behavior of $\mu$ has 
led \citep{Golden:PRL-3935,Murphy:JMP:063506}
%~\citep{Golden:PRL-3935,Murphy:JMP:063506}
to a detailed description
of these critical transitions in $\sigma^*,$ which is 
% directly 
analogous to
the Lee--Yang--Ruelle--Baker
description \citep{Baker-1990,Golden:PRL-3935} of the 
Ising model phase transition in the magnetization
$M.$ Moreover, using this gap behavior, all of the classical critical
exponent scaling relations
%~\citep{Bergman:SSP-147}
were recovered~\citep{Murphy:JMP:063506,Golden:PRL-3935} without heuristic
scaling forms \citep{Efros:PSS-475} but instead by using 
the \emph{rigorous} integral representation for
$\sigma^*$ involving $\mu$. 




This spectral behavior emerges in all the systems mentioned above,
% variety of composite systems, 
such as the brine microstructure of sea ice~\citep{Golden:S-2238,Golden:GRL:L16501,Golden:NAMS:2009} as shown in \figref{fig:brine_SM_E}, melt
ponds on the surface of Arctic sea ice~\citep{Hohenegger:TC-2012} as shown in \figref{fig:melt_ponds_SM_E}, and
the sea ice pack itself \citep{Murphy:PRL:2017} in Figure \ref{fig:pack_ice_SM_E}.
% (as shown in \figref{spectral_measure}a--c).
% and porous human bone~\citep{Kabel:B:1999}.
%
This also gives rise to critical behavior of the electric field, as shown in 
\figref{fig:brine_SM_E} for 2D cross sections of 3D brine microstructure, with $E_0$
taken to be in the vertical direction. Disconnected and weakly connected
examples of brine microstructure have small values of the electric field, while
strongly connected brine microstructures are characterized by a substantial 
increase the electric field's strength with "hot spots" forming in geometric bottlenecks. 
A similar behavior is exhibited by the 
temperature gradient $\nabla T$ associated with the Stieltjes integral for the effective
horizontal thermal conductivity $\kappa^*$, as shown for melt ponds atop Arctic sea ice in
\figref{fig:melt_ponds_SM_E} and for Arctic pack ice in \figref{fig:pack_ice_SM_E}.

\subsection{Generalization to rank deficient setting}
%
In the periodic setting, for example, the matrix Laplacian is singular so the matrix representation of 
%the operator 
$(-\Delta)^{-1}$ in $\Gamma$ is not defined. We now extend the mathematical framework developed in \citep{Murphy:2015:CMS:13:4:825}
to this setting. To make the connection to the abstract Hilbert space \citep{Golden:CMP-473} and full rank matrix \citep{Murphy:2015:CMS:13:4:825} settings, we first give relevant details for these cases. Equation \eqref{eq:Electric_Field_two_component} for the abstract Hilbert space setting follows by applying the operator $-\nabla(-\Delta)^{-1}$ to the formula $\nabla\cdot\vecD=0$, yielding $\Gamma\vecD=0$. Equation \eqref{eq:Electric_Field_two_component} then follows by using $\Gamma E_f=E_f$ and $\Gamma E_0=0$ \citep{Murphy:2015:CMS:13:4:825}, since $E_f$ is in the range of $\Gamma$ and $E_0$ is constant \citep{Murphy:ADSRF-2019,Murphy:ADSTPF-2017,Murphy:2015:CMS:13:4:825}. The matrix form of $\nabla\cdot D=0$ is $-\nabla^T D=0$, where $\nabla$ now represents the finite difference matrix representation of the gradient operator and $-\nabla^T$ is the finite difference representation of the divergence operator, with \emph{negative} matrix Laplacian given by $\nabla^T\nabla$ \citep{Murphy:2015:CMS:13:4:825}. As before, in \citep{Murphy:2015:CMS:13:4:825} we applied the matrix $\nabla(\nabla^T\nabla)^{-1}$ to the formula $-\nabla^T D=0$, yielding $\Gamma\vecD=0$, where $\Gamma=\nabla(\nabla^T\nabla)^{-1}\nabla^T$, and equation \eqref{eq:Electric_Field_two_component} follows the same way as before.





Now consider the singular value decomposition of the matrix gradient \citep{Murphy:ADSRF-2019} of size $m\times n$, say, $\nabla=U\Sigma V^T$. Here $U$ is a $m\times n$ matrix satisfying $U^T U=I_n$, $\Sigma$ is a $n\times n$ diagonal matrix with diagonal entries consisting of the \emph{singular values} of $\nabla$, and $V$ is a $n\times n$ orthogonal matrix satisfying $V^T V=V V^T=I_n$, where $I_n$ is the identity matrix of size $n$. When the matrix gradient is full rank it has $n$ strictly positive singular values, so $\Sigma$ is an invertible matrix and the matrix representation of $\Gamma$ is given by $\Gamma=UU^T$. On the other hand, when the matrix gradient is singular we have 
$\Sigma = \text{diag}(\Sigma_1,0,\ldots,0)$, where the diagonal matrix $\Sigma_1$ contains the $n_1$ strictly positive singular values of $\Sigma$ and the rest of the singular values have value 0. Denoting $U_1$ and $V_1$ to be the columns of $U$ and $V$ corresponding to the diagonal entries of $\Sigma_1$, we have $\nabla=U_1\Sigma_1 V_1^T$, where $\Sigma_1$ is invertible and $U_1^T U_1=V_1^T V_1=I_{n_1}$. This enables us to write $-\nabla^T D=0$ as $-V_1\Sigma_1 U_1^T D=0$, hence $U_1^T D=0$ and $U_1U_1^T D=0$. Noting that the columns of $U_1$ span the range of the matrix gradient $\nabla$, the matrix $U_1U_1^T$ is a projection onto the range of $\nabla$ \citep{Murphy:ADSRF-2019}. Defining $\Gamma=U_1U_1^T$, equation \eqref{eq:Electric_Field_two_component} follows the same way as before. This clearly generalizes the full rank setting. More details will be published elsewhere.

\section{Analytic continuation for polycrystalline media}
%
Sea ice is a composite material with polycrystalline 
% and porous brine microstructures
microstructure
on the millimeter to centimeter scale. When sea water freezes under
turbulent conditions, granular sea ice forms, which has small 
crystals with isotropic orientation angles. Columnar sea ice forms 
in quiescent conditions, with large crystals more strongly
oriented in the vertical direction. Examples of granular and 
columnar sea ice polycrystalline microgeometries are displayed 
in Fig. \ref{polycrystalline} (a).


Our analysis 
of the
% electromagnetic 
transport properties of random, uniaxial
polycrystalline media \citep{Barabash:JPCM:10323}
% has demonstrated that 
in \citep{Gully:PRSA:2015}, and a somewhat new formulation
presented below, shows
the underlying
mathematical framework is a direct analogue of that for two-phase
random media discussed in Sec. \ref{sec:two_component}.
For simplicity, we discuss electrical permittivity $\epsilon$, 
keeping in mind the broader applicability to thermal conductivity $\kappa$, 
electric conductivity $\sigma$, etc. Polycrystalline materials, 
are composed of many crystallites (single crystals 
of varying size, shape, and orientation) that can have different local
conductivities along different crystal axes.  In contrast to 
equation \eqref{eq:locally_iso_sigma}, the local
permittivity matrix of such media is given by
\citep{MILTON:2002:TC,Barabash:JPCM:10323}  
% $\beps(x,\omega)=R^{\,T}\,\Phi R$,  
%
\begin{align}\label{eq:polycrystal_parameters}
\beps(x,\omega)=R^{\,T}\,\Phi R\,,
\quad
\Phi = \text{diag}(\epsilon_1,\ldots,\epsilon_d)\,,
\end{align}
%
where $R(\vecx,\omega)$ is a random rotation matrix satisfying
$R^{\,T}=R^{\,-1}$. For example, for $d=2$ we have 
%
\begin{align}\label{eq:polycrystal_parameters_2D}
\beps=R^{\,T}
\left[
\begin{array}{cc}
	\epsilon_1& 0\\
	0 & \epsilon_2\\
\end{array}
\right]
R,
%     \quad
%     \brho=R^{\,T}
%   \left[
%   \begin{array}{cc}
	%     1/\epsilon_1& 0\\
	%     0 & 1/\epsilon_2\\
	%     \end{array}
%     \right]
%     R,
\qquad
R=
\left[
\begin{array}{rr}
	\cos\theta& -\sin\theta\\
	\sin\theta & \cos\theta\\
\end{array}
\right],
\end{align}
%
where $\theta=\theta(\vecx,\omega)$ is the orientation angle, measured from the
direction $\vece_1$, of the crystallite which has an interior 
containing $\vecx\in\mathbb{R}^d$ for $\omega\in\Omega$. 
% Here, $\vece_j$, $j=1,\ldots,d$, are standard basis vectors with components $(\vece_j)_k=\delta_{jk}$. 
% and $\delta_{jk}$ is the Kronecker delta. 
In higher
dimensions, $d\geq3$, the rotation matrix $R$ is a composition of
``basic'' rotation matrices 
$R_i$, e.g. $R=\prod_{j=1}^dR_j$, where the matrix $R_j(\vecx,\omega)$
rotates vectors in $\mathbb{R}^d$ by an angle
$\theta_j=\theta_j(\vecx,\omega)$ about the $\vece_j$ axis. For example, in three
dimensions 
%
\begin{align}\label{eq:polycrystal_parameters_3D}
R_1=
\left[
\begin{array}{ccc}
	1     &0     &    0\\
	0     &\cos\theta_1 & -\sin\theta_1\\
	0     &\sin\theta_1 & \cos\theta_1\\
\end{array}
\right],
\quad
R_2=
\left[
\begin{array}{ccc}
	\cos\theta_2  & 0     &\sin\theta_2 \\
	0      & 1     & 0\\
	-\sin\theta_2 & 0     &\cos\theta_2\\
\end{array}
\right],
\quad
R_3=
\left[
\begin{array}{ccc}
	\cos\theta_3 & -\sin\theta_3 & 0\\
	\sin\theta_3 & \cos\theta_3  & 0\\
	0     & 0      & 1\\
\end{array}
\right].
\end{align}
%



In the case of \emph{uniaxial} polycrystalline media, the local permittivity along one of the crystal
axes has the value $\epsilon_1$, while the permittivity along
all the other crystal axes has the value $\epsilon_2$, so
$\Phi=\text{diag}(\epsilon_1,\epsilon_2)$ for 2D (which is the general setting for 2D) 
and $\Phi=\text{diag}(\epsilon_1,\epsilon_2,\epsilon_2)$ for 3D. 
Equation~\eqref{eq:polycrystal_parameters} can be written in a more
suggestive form in terms of the matrix $C=\text{diag}(1,0,\ldots,0)$  
%
\begin{align}\label{eq:sigma_uniaxial}
\beps(\vecx,\omega)=\epsilon_1X_1(\vecx,\omega)+\epsilon_2X_2(\vecx,\omega), \qquad
%   \brho(\vecx,\omega)=X_1(\vecx,\omega)/\epsilon_1+X_2(\vecx,\omega)/\epsilon_2.
%   \beps=\epsilon_1X_1+\epsilon_2X_2,
\end{align}
%
which is an analogue of equation \eqref{eq:sigma_two-phase}.
Here $X_1=R^{\,T}CR$ and $X_2=R^{\,T}(I-C)R$, where $I$ is the
identity matrix on $\mathbb{R}^d$. Since $R^{\,T}=R^{\,-1}$ and $C$
is a diagonal projection matrix satisfying $C^{\,2}=C$, it is
clear that the $X_i$, $i=1,2$, are mutually orthogonal projection
matrices satisfying
%$X_jX_k=X_j\delta_{jk}$ and $X_1+X_2=I$.
%
\begin{align}\label{eq:Projection_Matrices}
X_j^{\,T}=X_j, \quad X_jX_k=X_j\delta_{jk}, \quad X_1+X_2=I,
\end{align}
%
which are also properties of the characteristic functions $\chi_j$ in Sec. \ref{sec:two_component}.


Equations \eqref{eq:Maxwells_Equations} and \eqref{eq:E_variational} are also satisfied in 
this polycrystalline setting \citep{Golden:CMP-473}. Similar to the derivation of 
equation \eqref{eq:Electric_Field_two_component} in Sec. \ref{sec:two_component}, a resolvent 
representation for $X_1E$ follows by applying the operator $-\nabla(-\Delta)^{-1}$ to the 
formula $\nabla\cdot\vecD=0$, yielding $\Gamma\vecD=0$. Then, using $\Gamma E_f=E_f$ and 
$\Gamma E_0=0$ \citep{Murphy:2015:CMS:13:4:825} yields
the following analogue of equation \eqref{eq:Electric_Field_two_component} 
%
\begin{align}\label{eq:Electric_Field_unaxial}
X_1\vecE=s(sI-G)^{-1}X_1\vece_k\,,
%\chi_1\vecE=s(sI-\chi_1\Gamma\chi_1)^{-1}\chi_1e_k\,,
%
% \chi_1\vecE=(sI-G)^{-1}\xi\,,
\quad
G=X_1\Gamma X_1,
\end{align}
%
%\vspace{-0.5ex}
yielding the integral representation in equation \eqref{eq:F(s)} for 
$F(s) = \langle [(sI - G)^{-1} X_1 \vece_k] \cdot \vece_k \rangle$.
As in the two component setting, a critical feature of 
equation \eqref{eq:F(s)} is that
the component parameters 
in $s$ are separated from the geometrical information in $\mu$.
% (Extensions to multicomponent media involve several complex variables
% \citep{Golden:JSP-655,Golden:JMPS-333,Milton:CPAM-647,Milton:CMP-281,DellAntonio:JMP-2688}.)
Information about the geometry enters
through the moments in equation \eqref{eq:Moments_Fs} with $G$ given in 
\eqref{eq:Electric_Field_unaxial} and $\chi_1$ replaced by $X_1$.
The mass $\mu_0$ of the measure $\mu_{jk}$ is given by    
%
\begin{align}\label{eq:Mass_mu}
\mu_{jk}^0=\langle X_1\vece_j\cdot\vece_k\rangle, \quad
\mu_{kk}^0=\langle |X_1\vece_k|^2\rangle,
\end{align}
%
where the second equality follows from the fact that $X_1$ is a
real-symmetric projection matrix.  
The statistical average $\langle |X_1\vece_k|^2\rangle$ in~\eqref{eq:Mass_mu}
can be thought of as the ``mean orientation,'' or as the percentage of
crystallites oriented in the $k^{\text{th}}$ direction.
% Since the \emph{projection} 
% matrix $X_1$ is bounded by one in operator norm, the
% Cauchy--Schwartz inequality and~\eqref{eq:Mass_mu} imply that
% $0\leq\mu_{kk}^0\leq1$. 
For example, in the case of  
two-dimensional polycrystalline media, $d=2$, equation
\eqref{eq:polycrystal_parameters_2D} implies that 
%$\mu_{11}^0=\langle\cos^2\theta\rangle$, $\mu_{22}^0=\langle\sin^2\theta\rangle$, and
%$\mu_{12}^0=\langle\sin\theta\cos\theta\rangle$.
%
\begin{align}
\mu_{11}^0=\langle\cos^2\theta\rangle, \quad
\mu_{22}^0=\langle\sin^2\theta\rangle, \quad
\mu_{12}^0=\langle\sin\theta\cos\theta\rangle.
\end{align}
%
Generalizing equation~\eqref{eq:polycrystal_parameters_3D}, with
$R=\prod_{j=1}^dR_j$, to dimensions $d\geq3$ shows that $\mu_{jk}^0$ is a
linear combination of averages of the 
form $\langle\prod_i\cos^{n_i}\theta_i\sin^{m_i}\theta_i\rangle$, where $n_i,m_i=0,1,2,\ldots$.



Given partial information on the 
microgeometry incorporated into the moments $\mu_n$, the integral representation \eqref{eq:F(s)} for this polycrystalline setting
yields rigorous \emph{forward bounds} for the 
effective parameters of composites \citep{Gully:PRSA:2015,MILTON:2002:TC}, as shown
in \figref{polycrystalline}d below.
The integral representation can also be used to obtain \emph{inverse
bounds} that give estimates for the structural parameters of a composite derived from  
the electromagnetic response
of a sample. One such structural characteristic is 
the average crystallite orientation \citep{Gully:PRSA:2015,MILTON:2002:TC}; see 
Section \ref{inverse} for more details.





\begin{figure}%[h!] 
\centering
\includegraphics[width=\textwidth]{polycrystal_E_SM_Bounds.png}
\caption{\emph{Spectral analysis of polycrystalline media.} 
	(a) Cross sections of polycrystalline microstructure for granular and columnar sea ice (J.-L. Tison). (b) Discrete checkerboard polycrystal microstructure with isotropic crystallite 
	orientations within the horizontal plane, with small (top) and large (bottom) crystallite size. Cool and warm colors correspond to low and high displacement field values. (c) The spectral function, a histogram representation of the spectral 
	measure $\mu(\lambda)$ shown along with it's theoretical prediction for such 
	isotropic media \citep{MILTON:2002:TC}. (d) An example value of the complex
	effective permittivity of isotropic polycrystalline media captured by first and second order bounds \citep{Gully:PRSA:2015}. 
 %%%   \textcolor{red}{From Daniel: is there a reason the bounds in (d) are so much bigger than the bounds in figure 7? I do not see this plot in \citep{Gully:PRSA:2015}}
}
\label{polycrystalline}      
\end{figure}



\subsection{Spectral measure computations for uniaxial polycrystalline materials} 
\label{polycrystal_discrete}
%
Computing the spectral measure $\mu$ for a given polycrystalline microgeometry first
involves discretizing the composite into a square lattice with vertex values in the range 
$[0,2\pi]$ corresponding to the crystallite orientation angles at each vertex location.
On this square lattice
the action of the differential operators $\bnabla$ and $\bnabla\cdot$ are 
defined in terms of 
forward and backward difference operators
\citep{Golden:CMP-467, Murphy:2015:CMS:13:4:825}.
Then the key operator 
$X_1\Gamma X_1$, 
which depends on the geometry of the 
network via $X_1$, becomes a 
%self adjoint matrix. 
%nbmurphy
real-symmetric matrix $M$.
%
Here $\Gamma$ is as in Sec. \ref{sec:two_component_discrete}
and $X_1$ is a banded (random) projection matrix that 
determines the geometry of the polycrystalline medium. 
% In \figref{fig:Gamma_Matrix_2D_3D} we visually display the matrix entries of $\Gamma$ with periodic boundary conditions. 
% for dimension $d=2$ and $d=3$. 
% The clear block matrix structure determines the properties of the components $\mu_{jk}$ of the matrix valued spectral measure $\bmu$ associated with anisotripic random media \citep{Murphy:2015:CMS:13:4:825}.
In this setting, the integral and $X_1\vecE$ in equations \eqref{eq:F(s)} and \eqref{eq:Electric_Field_two_component} 
have explicit representations in terms of the eigenvalues $\lambda_i$ and eigenvectors $\vecw_i$ of $M$ \citep{Murphy:2015:CMS:13:4:825} 
given by equation \eqref{eq:Discrete_E_F(s)}, and similarly the spectral measure is given by equation \eqref{eq:Stieltjes_F_Discrete_two-component}, with $\chi_1$ replaced by $X_1$. 
%

% We have recently developed numerical methods 
% analogous to that in Section \ref{Math_models}, 
% enabling the direct computation 
% of spectral measures, 
% bulk transport coefficients, local fields, 
% eigenvalue statistics, and eigenvector and field localization properties,
% as displayed in \figref{polycrystalline}. 
% We have recently computed the spectral measure $\mu$ and the effective thermal and electrical conductivities, $\kappa^*$ and $\epsilon^*$, via \eqref{eq:Stieltjes_F_Discrete_two-component} for 
% % polycrystalline 
% microscructure corresponding to thin 2D  polycrystaline films with random checkerboard microgeometry and isotropic orientations. 
% We observe that $\kappa^*$ and $\epsilon^*$ increase with increasing grain size, and also the formation of significant energy barriers along grain boundaries, which increase in intensity with increasing grain size, which is consistent with other
% experimental \citep{Liu:2014:JPhysChemC:118:24797,Kuhn:2018:SciRep:8:14868,Limbu:2017:Carbon:117:367,Ma:2017:NatureComm:8:14486,Isacsson:2016:2DMat:4:1:012002,Ekimov:2011:InorgMat:44:224}, theoretical \citep{Luo:2009:ChinPhysLett:26:11,Dong:2014:SciRep:4:7037,Isacsson:2016:2DMat:4:1:012002}, and numerical \citep{Limbu:2017:Carbon:117:367,Isacsson:2016:2DMat:4:1:012002,Mortazavi:2015:SciRep:5:13228,Liu:2014:JPhysChemC:118:24797} methods and results. 








The following theorem provides a rigorous mathematical formulation of
integral representations for the effective parameters for finite
lattice approximations of random uniaxial polycrystalline media,
analogous to Theorem \ref{thm:Discrete_Spectral_Theorem_ACM_two-comp} above.    
This theorem formulated for the polycrystal problem, holds for both of the settings where the matrix gradient
is full rank or rank deficient. The proof of the theorem will be published elsewhere. 
% The
% theorem and proof are formulated in terms of the random matrix
% $M=X_1\Gamma X_1$. The formulations involving the matrices $M_2=X_2\Gamma X_2$
% and $K_i=X_i\Upsilon X_i$, $i=1,2$, are analogous.  

% \vspace{0.15in}
% 
\begin{theorem} \label{thm:Discrete_Spectral_Theorem_ACM_polycrystal}
For each $\omega\in\Omega$, let $M(\omega)=W(\omega)\Lambda(\omega)\,W(\omega)$ be the eigenvalue
decomposition of the real-symmetric matrix
$M(\omega)=X_1(\omega)\,\Gamma\,X_1(\omega)$. Here, the columns of the matrix $W(\omega)$
consist of the orthonormal eigenvectors $\vecwt_i(\omega)$, $i=1,\ldots,N$,
of $M(\omega)$ and the diagonal matrix $\Lambda(\omega)={\rm diag}(\lambda_1(\omega),\ldots,\lambda_N(\omega))$
involves its eigenvalues $\lambda_i(\omega)$. 
Denote $Q_i=\vecwt_i\,\vecwt_i^{\;T}$ the projection matrix onto the eigenspace
spanned by $\vecwt_i$.
The electric field $\vecE(\omega)$ satisfies $\vecE(\omega)=\vecE_0+\vecE_f(\omega)$, with
$\vecE_0=\langle\vecE(\omega)\rangle$ and $\Gamma\vecE(\omega)=\vecE_f(\omega)$, and the
effective complex permittivity tensor $\beps^*$ has components
$\epsilon^*_{jk}$,
%   $\epsilon_{jk}^*=\epsilon_2\,m_{jk}(h)$, 
$j,k=1,\ldots,d$,  which satisfy       
%
\begin{align}\label{eq:Stieltjes_F_Discrete_polycrystal}
	&\epsilon^*_{jk}=\epsilon_2(\delta_{jk}-F_{jk}(s)), 
	&&F_{jk}(s)=\int_0^1\frac{\d\mu_{jk}(\lambda)}{s-\lambda}\,, 
	&&\d\mu_{jk}(\lambda)=\sum_{i=1}^N\langle \delta_{\lambda_i}(\d\lambda)\;X_1\,Q_i\hat{e}_j\cdot\hat{e}_k\rangle\,.  
\end{align}
% 
\end{theorem}


To numerically compute $\mu$, a non-standard generalization of the spectral theorem for
matrices is required due to the projective nature of the matrices $X_1$ and $\Gamma$
\citep{Murphy:2015:CMS:13:4:825}. In particular, a
\emph{projection method} analogous to that in \citep{Murphy:2015:CMS:13:4:825} shows the spectral measure $\mu$ in
\eqref{eq:Stieltjes_F_Discrete_polycrystal} depends only on the eigenvalues and eigenvectors of the 
upper left $N_1\times N_1$ block of the matrix $R\Gamma R^T$, where $N_1=N/d$. These
submatrices are smaller by a factor of $d$, which improves the efficiency and 
numerical computations of $\mu$ by a factor of $d^3$.
In \figref{polycrystalline}, computations of the 
displacement field $\vecD$ are displayed for 2D
polycrystalline media for small and large crystal sizes alongside
cross sections of polycrystalline microstructure for granular
and columnar sea ice. When the effective permittivity tensor 
$\beps^*$ is diagonal, such as the setting of isotropically oriented
crystallites, the spectral measure for an infinite system is known
in closed form \citep{MILTON:2002:TC} to be 
$d\mu(\lambda)=(\sqrt{(1-\lambda)/\lambda})(\d\lambda/\pi)$, as
shown in \figref{polycrystalline} (c). This
measure has a singularity at $\lambda=0$, which indicates that the
material is electrically conductive on macroscopic length scales
\citep{Murphy:2015:CMS:13:4:825,Murphy:JMP:063506}. When 
the polycrystalline material is isotropic, 
both the mass and first moment of the measure $\mu$ are
known, which enables two nested bounds for $\epsilon$ to be computed 
\citep{Gully:PRSA:2015}, as shown in \figref{polycrystalline} (d).
%\textbf{Should we say anything about the transitional behavior as a function of crystallite size?}


%  \input{inverse}




\section{Inverse homogenization: Recovery of information 
about the microstructure of composites}
\label{inverse}

Developed originally for the effective complex permittivity $\e^*$, the integral representation (\ref{eq:F(s)}) 
%%
%\begin{equation} \label{eq:F(s)}
%    F(s)=1 - \frac{\e^*}{\e_2} =\int_0^1 {\frac{d\mu(\lambda)}{s-\lambda}}\,,
%   \qquad s = \frac{1}{1-\e_1/\e_2}\,,
%\end{equation}
%%
yields rigorous \emph{forward bounds} for the effective permittivity $\e^*$ of two-component composites formed of materials with permittivity $\e_1$ and $\e_2$, given partial information on the 
microgeometry via the moments $\mu_n$ 
\citep{Bergman:PRL-1285,Milton:APL-300,Bergman:AP-78,Golden:CMP-473}.
One can also use the integral representation to recover information about the structure of the composite material $-$ this is the  problem of \textbf{inverse homogenization} 
\citep{Cherkaev:IP-1203}. For inverse homogenization, it is important that the representation  (\ref{eq:F(s)}) separates information about the properties of the phases contained in the parameter  $s$  from information about the microgeometry contained in the measure $\mu$ and its moments $\mu_n=\langle G^{\,n} \chi_1 \vece_k \cdot \chi_1\vece_k \rangle$ in  (\ref{eq:Moments_Fs}) via higher-order correlation functions of the geometry function $\chi_1$. 


%A spectral measure is the key mathematical object appearing in integral representations for effective transport and mechanical properties of two phase composites. 

In principle, the spectral measure $\mu$ and its moments $\mu_n$ contain all the geometrical information about the composite. For example, the mass $\mu_0$ is the 
volume fraction $\phi$ of the first component in the composite,
\begin{equation}                                   \label{volfrac}
	\mu_0 = \int^1_0  d \mu (z) = 
	\langle \chi_1 \rangle = \phi,
\end{equation}
and the fraction of the second phase
is $1-\phi$. Connectivity information is also embedded in the spectral measure. 

%$$
%\mu_n = 
%\int^1_0 z^n d \mu (z)  =
%(-1)^n \langle \chi_1 [( \Gamma \chi_1 )^n \vece_k ] \cdot \vece_k \rangle.
%$$
%%
%\begin{align}\label{eq:Moments_Fs}
%\mu_n = 
%\int^1_0 \lambda^n d \mu (\lambda)  =
%%\langle [\chi_1\Gamma \chi_1]^{\,n} \vece_k \cdot \vece_k \rangle,
%\langle G^{\,n} \chi_1 \vece_k \cdot \chi_1\vece_k \rangle.
%%\quad
%%G=\chi_1 \Gamma \chi_1. 
%\end{align}
%%

The basis for inverse homogenization is provided by the uniqueness theorem \citep{Cherkaev:IP-1203}, which formulates the conditions under which the measure $\mu$ in the representation  (\ref{eq:F(s)}) can be uniquely reconstructed from measured data. 
For instance, complex permittivity data measured for a range of frequencies of the applied electromagnetic field is sufficient to uniquely recover the measure $\mu$ in (\ref{eq:F(s)}) \citep{Cherkaev:IP-1203}. Such data is also sufficient for the unique reconstruction of the moments $\mu_n$  \citep{Cherkaev:IP-065008}, provided the permittivity of one of the phases is frequency dependent. 
Two major approaches to inverse homogenization are (i) \emph{reconstruction of the measure $\mu$ }
\citep{Cherkaev:IP-1203,Cherkaev:IP-065008,Day:JPCM-96,Zhang:JCP-5390,Bonifasi-Lista:PMB-3063,Bonifasi:2009b,Cherkaev:JBiomech-345,Day:1999,Day:2000,Golden:JBM:337,Cherkaev:Royal:2020}  
and then calculating its moments, and (ii) \emph{inverse bounds} for the structural parameters, such as the volume fractions of the components    
\citep{McPhedran:APA:19,McPhedran:MRSP:1990:195,Cherkaev:1996,Cherkaev:WRM-437,Cherkaev:IP-1203,Cherkaev:IP-065008},
orientation of the crystals \citep{Gully:PRSA:2015}, or connectedness \citep{Orum:PRSLA:2012} of one phase. 



When only a few data points are available, though the uniqueness theorem \citep{Cherkaev:IP-1203} is not immediately applicable, one can outline a set of measures consistent with the measurements, 
\begin{equation}
	{\cal M} = \{ \mu: F_{\mu}(s) = 
	1 -\epsilon^\ast / \epsilon_2\} . 
	\label{setM}
\end{equation} 
Then determine an interval confining the first moment of the measure  $\mu$ providing, for instance,  an interval of uncertainty for the volume fraction of one material. For several data points corresponding to the same structure of the composite, e.g., measurements at distinct frequencies, the bounds for the volume fraction are given by an intersection of all admissible intervals \citep{Cherkaev:1996,Cherkaev:WRM-437,Tripp:1998}. 
When the requirements for the measurements needed to uniquely reconstruct the spectral measure $\mu$ established by the uniqueness theorem are satisfied, the set ${\cal M}$ is reduced to one point. 
But the map from the set of measures to the set of microstructures is not unique, and  
a variety of microgeometries generate the same response under the applied 
field. Different microgeometries corresponding to the same sequence of moments $\mu_0, \mu_1, ... $ are the $S-$equivalent structures \citep{Cherkaev:IP-1203} that are not distinguishable by homogenized measurements. 

%  

An equivalent representation for the function $F(s) $ in \eqref{eq:F(s)} using a logarithmic potential of the measure $\mu$ in the complex $s-$plane is  \citep{Cherkaev:IP-1203}:
\begin{equation}                               \label{ln}
	F (  s  )   =   
	\, \frac{d}{ds} \,  % \frac{\D} {\D s }\, 
	\int \ln\,  | s - z | \, \,   d \mu(z)  , \qquad 
	   d / ds  =   (\D / \D x - i~  \D/\D y)  \;, 
	\qquad s = \frac{1}{1-\e_1/\e_2}\,.
\end{equation}
%\begin{equation}                               \label{ln}
%F (  s  )   =
%   \,  \frac{\D \Phi} {\D s }\,  , \qquad
% \Phi(s) = \int_0^1 \ln\,  | s - \lambda | \, \,   d \mu(\lambda) \;,
%\end{equation}
%where $\frac{\D}{\D  s} =   \frac{\D}{\D x} - i \frac{\D}{\D y}\;$ 
The solution to the inverse problem of recovering the measure $\mu$ is constructed by solving the minimization problem: 
\begin{equation}
	\mbox{min}_{\mu} \; || A \mu - F ||^2 \,, \qquad   F(s)= 1-\e^*(s)/\e_2\;.
\end{equation}
Here $A$ %$A\mu = \frac{\D}{\D s}\int_0^1 \ln\, | s - \lambda | \, \,d \mu(\lambda) $, 
is the integral operator in  (\ref{ln}) or in (\ref{eq:F(s)}), the norm is the $L^2-$norm, $F$ is the given function of the measured data, and $s \in \mathcal{C}$, where $\mathcal{C}$  is a curve in the complex plane corresponding to the range of frequencies of the applied field.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The solution of the minimization problem does not depend continuously on the data. Unboundedness of the operator $A^{-1}$  leads to arbitrarily large variations in the solution, and the problem requires regularization to design a stable numerical algorithm \citep{Cherkaev:IP-1203}. 
Regularized inversion schemes and stable reconstruction algorithms  to recover $\mu$ and its moments from data on the effective complex permittivity
were developed in \citep{Cherkaev:IP-1203,Cherkaev:2004:331,Cherkaev:IP-065008,Bonifasi-Lista:PMB-3063,Cherkaev:JBiomech-345}. They are
based on $L^2$, total variation ($TV$), non-negativity constraints, and constrained Pad\'e approximations of the measure $\mu$ \citep{Zhang:JCP-5390}.
%
In applications to imaging of bone structure, spectral measures $\mu$ are computed with regularization algorithms based on $L^2$ constrained minimization, from electromagnetic \citep{Bonifasi-Lista:PMB-3063,Cherkaev:JBiomech-345,Golden:JBM:337} and viscoelastic \citep{Bonifasi:2008,Bonifasi:2009b,Cherkaev:JBiomech-345} data. 
This allows one to separate samples of healthy and osteoporotic bone by distinguishing the relative fractions of bone tissue and marrow in the different microstructures, as well as the connectivity of the trabecular architectures.





%%************    Generalizations: Viscoelastic   *******
%Viscoelasticity: uniaxial polycrystalline material 

The first application of Stieltjes representations to 
the elastic properties of two-phase composites can be found in
\citep{Kantor:1982:JMPS,Kantor:1984:JMPS}.
Now consider hydrostatic and deviatoric projections $\Lambda_h$ and $\Lambda_s$ onto
the orthogonal subspaces of the second order tensors which are
proportional to the identity tensor and trace-free. 
Then the 
Stieljtes integral representation
was generalized in \citep{Cherkaev:JBiomech-345} to the effective viscoelastic modulus 
% for arbitrary loadings 
and to two-dimensional viscoelastic polycrystalline materials \citep{Cherkaev:AAM:2019} under the assumption that the constituents have the same elastic
bulk and different (elastic and viscoelastic) shear moduli. This representation was also 
used in inverse homogenization \citep{Bonifasi:2008,Cherkaev:JBiomech-345,Cherkaev:Royal:2020} for successful recovery of the volume fractions of the phases in a composite from known viscoelastic shear moduli. 

Other approaches to the volume fraction bounds include 
\citep{Engstrom:2005,Milton:2012universal,Thaler:2014},
based on estimates for higher order moments and on variational bounds, as well as direct inversion of known formulas or mixing rules \citep{Bergman:SSP-147,Levi:2013} for effective properties of composites with specific structure. However, an advantage of the methods discussed here is their applicability without a priori assumption about the microgeometry.
% *************
%

{\bf Spectral coupling of various properties of composites.}
%
An important application of inverse homogenization is for indirectly evaluating material properties through cross-coupling \citep{Cherkaev:IP-1203,MILTON:2002:TC}. Different properties of composites are coupled through their microgeometry; this relationship can be used for estimating properties from available data that are otherwise difficult to measure. The conventional approaches are based on empirical and semi-empirical relations, such as, for instance, Kozeny-Carman or Katz-Tompson. These relations estimate the permeability of a porous material characterizing the microstructure by a ''formation factor'' $F$, which relates properties of one phase in the composite to the effective properties 
of the material \citep{Sahimi-95,Torquato:RHM-02,Wong:PRB-6606,Wong:PT:1988}.
% $F = \s/\s^*$, where  $\s$ is the conductivity of a fluid filling in the porous material, and $\s^*$ is the effective conductivity.  

In the \emph{spectral coupling} method \citep{Cherkaev:IP-1203} based on the Stieltjes representation (\ref{eq:F(s)}), the spectral measure $\mu$ associated with the geometric structural function, couples various properties of the same material. Spectral coupling \citep{Cherkaev:IP-1203,Cherkaev:2004:331,Cherkaev-Zhang:2003,Cherkaev:JBiomech-345} for two-component composites allows us to recover various transport properties of sea ice from the spectral measures computed using other measured properties. In particular, this approach can utilize effective complex permittivity data (recovered from radar measurements) to calculate the thermal conductivity \citep{Cherkaev-Zhang:2003} and hydraulic conductivity of polycrystalline sea ice, which are normally difficult to measure over large scales.
Spectral coupling was also extended to evaluating viscoelastic properties of two-component composites in \citep{Cherkaev:JBiomech-345}, for applications to characterizing bone properties and microarchitecture. 

Inverse homogenization for recovering microstructural parameters from effective property measurements is useful in multiple forms of non-destructive testing of materials, specifically in applications such as medical imaging as well as Synthetic Aperture Radar (SAR) remote sensing for assessing the structure and transport properties of sea ice.

\subsection{Bounds for the moments of the spectral measure}
%%%%%%%%%%%%%%%%%%%%%%%


The second approach to the inverse homogenization problem is calculating \emph{inverse bounds} for the structural parameters, such as the volume fraction of each of the components    
\citep{McPhedran:APA:19,McPhedran:MRSP:1990:195,Cherkaev:1996,Cherkaev:WRM-437,Cherkaev:IP-1203}, 
orientation of the crystals \citep{Gully:PRSA:2015}, or connectedness \citep{Orum:PRSLA:2012} of, say, the conducting phase. 
An analytical approach to estimating the volume fractions of materials in a composite \citep{Cherkaev:1996,Cherkaev:WRM-437,Tripp:1998} gives explicit analytic formulas for the first order inverse bounds on the volume fractions of the constituents in a general composite and second order inverse bounds on the fractions of the phases in an isotropic composite \citep{Cherkaev:WRM-437}.


The inverse bounds are derived using analyticity of the effective complex permittivity of the composite. The first order bounds  $p_l^{(1)}$ and $p_u^{(1)}$ for the volume fraction $\phi$ give the lower and upper bounds for the zeroth moment $\mu_0$ of the measure $\mu$ or its mass in (\ref{volfrac}) \citep{Cherkaev:1996,Cherkaev:WRM-437}:  
\begin{equation} 
	p_l^{(1)}  \le \phi \le p_u^{(1)}, 
	\qquad  
	p_l^{(1)}  
	= |f|^2 \, \frac{ Im\,(\bar{s}) }{Im\,(f)}  , 
\qquad  
p_u^{(1)}  
= 1 - \frac{ |g|^2 \, Im\,(\bar{t}) }{Im\,(g)}. 
\label{inv1b}
\end{equation}
Here, $t=1-s$, $f$ is the known value of $F(s)$, and $g$ is the known value of $G(t)=1-\epsilon^*/\epsilon_1$.


%  Fig INVERSE  %
\begin{figure}
\centering
\includegraphics[width=\textwidth]{inverse_bounds_hq.png}
%\includegraphics[width=2.5in]%[width=\textwidth]
%{inv-f2.pdf} ~~ 
%&
%\includegraphics[width=1.3in]{inv-q-separation}
%&
%\includegraphics[width=2.45in]{inv-polycrystal-forw}
%&    a \hspace{1.8in} & b \hspace{2.2in} & c
\caption{\emph{Forward and Inverse bounds. }
	(a) Illustration of bounds on the volume fraction of one component
	in the mixture derived from first order anisotropic  bounds (left panel), 
	and from the second order isotropic  bounds (right panel) for the  
	effective permittivity \citep{Cherkaev:WRM-437}. 
	The small lens shaped domains each contain $\e^*$ of the anisotropic (left) 
	and isotropic (right) composites corresponding to the volume fractions of the first component  
	$p_l$ and $p_u$ which give the lower and upper bounds for  
	the fraction of the first material.  
	(b) Lower bounds on the separation parameter $q_{min}$ versus temperature \citep{Orum:PRSLA:2012} calculated using data on the effective complex permittivity. 
	The inverted data clearly indicate that as the ice warms, the separations of the brine inclusions decrease. Stars and squares indicate different sea ice slabs. 
	(c) Polycrystalline bounds for different brine volume fractions \citep{Gully:PRSA:2015} for the permittivity of sea ice (left) together 
	with the measured effective permittivity of sea ice in \citep{Arcone:JGR-14281}. Comparison between two-component and polycrystalline bounds (right) $-$ the forward polycrystalline bounds shown in blue are a dramatic improvement over the two-component forward elementary bounds and isotropic bounds shown in red.
}
\label{inv-bounds}
\end{figure}
% END  Fig INVERSE   

First and second order forward and inverse bounds are illustrated in Fig. \ref{inv-bounds}(a)  \citep{Cherkaev:WRM-437}
where first order bounds for the effective complex permittivity of all anisotropic composites that could be formed from two materials of permittivity $\e_1$ and  $\e_2$ are presented in the left panel, while the second order isotropic bounds are shown in the right panel.  
The small lens shaped domains each contain the anisotropic (left) and isotropic  (right)  
mixtures corresponding to the volume fractions $\phi$ of the first component equal to $p_l^{(q)}$ and $p_u^{(q)}$, $q = 1,2$. The points $p_l^{(q)}$ and $p_u^{(q)}$ give the lower and upper bounds for the volume fraction of the first material in the composite. Superscripts $q=1$ and $q=2$ indicate the first and second order bounds. 
%
For a set of data points $\epsilon^*(k)$, $k=1,...,N$ corresponding to composites with the same structure, 
the bounds for the fraction $\phi$ of the first phase in the material are given by an intersection of all admissible intervals 
$p_l^{(q)}(k) \leq \phi \leq p_u^{(q)}(k)$ \citep{Cherkaev:1996}: 
\begin{equation}
P_l^{(q)} = \max_k \, p_l^{(q)}(k)  \, \le \, \phi  
\, \le \, \min_k \, p_u^{(q)} (k) 
= P_u^{(q)}, \qquad q = 1,2.
\end{equation}
\noindent
Here, $  p_l^{(q)} (k)$  and $  p_u^{(q)} (k) $ are lower and upper  bounds  for the volume fraction  $\phi$ 
calculated using the effective complex permittivity $\epsilon^\ast(k)$, respectively; $q$ is the order of the bounds, with $q=1$ for a general mixture and 
$q=2$ for an isotropic composite. 

In \citep{Cherkaev:WRM-437}, this method was applied to estimating brine volume fraction in sea ice from two data sets of 4.75 GHz measurements of the complex permittivity $\e^*$ of sea ice \citep{Arcone:JGR-14281} at $-6^\circ$C and  $-11^\circ$C  with fractions of brine $\phi=0.036$ and  $\phi=0.0205$. Sea ice was considered as a composite of three components: pure ice, brine, and air. The effective complex permittivity of the mixture of ice and air was calculated with the Maxwell-Garnett formula. The first order bounds estimate the brine  volume fraction as 
$ 0.0213 \le \phi \le 0.0664 $ and 
$ 0.0119 \le \phi \le 0.0320 $, for the data set 1 and 2, respectively. 
The second order inverse bounds derived 
% from the measurements of the effective complex permittivity 
with the assumption of 2D isotropy in the horizontal plane give the following estimates for the brine volume fraction: 
$ 0.0333 \le \phi \le 0.0422 $ for the first data set with brine volume $\phi=0.036$,  
and $ 0.0189 \le \phi \le 0.0213 $ for the second data set with volume fraction of brine $\phi=0.0205$. The first order bounds are further extended to polycrystalline materials and allow estimation of the mean crystal orientation \citep{Gully:PRSA:2015}.



\subsection{Matrix particle forward and inverse bounds}

Another parameter important in characterizing the structure of a composite material consisting of inclusions within a host matrix is the separation between inclusions. 
Inclusion separation is an indicator of the connectedness of phases -- a key feature in critical behavior and phase transitions; the separation parameter may be used to estimate closeness to the percolation phase transition.

Composites with non-touching inclusions of one material embedded in a host matrix of a different material are called matrix particle composites. For a matrix particle composite with
separated inclusions, tighter bounds on the effective complex permittivity may be obtained. 
In \citep{Orum:PRSLA:2012}, sea-ice is considered as a matrix particle composite in which the brine phase is contained in separated, circular discs of radii $r_b$ randomly located on a horizontal plane, each surrounded by a ``corona'' of ice with outer radius $r_i$. Such a material  
is called a $q$-material, where $q= r_b/r_i$. The minimal separation of brine inclusions is $2(r_i-r_b)= 2 r_i (1-q)$. 
In this case, as it is shown in \citep{Bruno:PRSLA-353}, the support of $\mu$ in (\ref{eq:F(s)}) 
lies in an interval $[s_m,s_M]$, $0<s_m<s_M<1$  %\citep{Bruno:91}.  
such that   
$s_m = \frac{1}{2}(1-q^2 )$ and $s_M=\frac{1}{2}(1+ q^2)$. 
%
The further the separation of the inclusions, the smaller the  
interval $[s_m,s_M]$, and the tighter the bounds. 
Smaller $q$ values indicate well separated brine 
(and colder temperatures as in Fig. \ref{inv-bounds}), 
while $q=1$ corresponds to no restriction on the separation, 
with $s_m=0$ and  $s_M=1$.

Two parameters characterizing the structure of the sea ice composite are the volume fraction $\phi$ of the brine inclusions and a separation parameter $q$ that quantifies how close the inclusions are to each other.  
Using observed values of effective complex permittivity  
and inverting the forward matrix particle bounds, information about these two parameters is obtained in \citep{Orum:PRSLA:2012} by  
solving a reduced inverse spectral problem and bounding the volume fraction of the constituents.
Inverse bounds for inclusion separation are shown in Fig. \ref{inv-bounds} 
\citep{Orum:PRSLA:2012}, where the lower bound $q_{min}$ is displayed versus temperature of the sea ice slab. The inverted data clearly indicate that as the
ice warms, the separations of the brine inclusions decrease. 
It is remarkable that this important phenomenon is characterized from electromagnetic measurements through an inversion scheme. %Stars, slab 84-3; squares, slab 84-4.


\subsection{Extension to polycrystalline composites}

The method of inverse bounds \citep{Cherkaev:1996,Cherkaev:WRM-437,Tripp:1998} for structural parameters of a composite from measured effective properties was extended to polycrystalline materials in \citep{Gully:PRSA:2015}.  In the case of a uniaxial polycrystalline composite,  \cite{Gully:PRSA:2015} developed bounds for the mean orientation of crystals in the 
sea ice from measured values of ice permittivity. 
As columnar and granular microstructures
have different mean single crystal orientations \citep{Weeks:CM-82}, % \citep{Weeks:GSPSI-82}, 
this inverse approach 
is useful for determining ice type when using remote sensing techniques.

The structures of different types of ice formed under different 
environmental conditions vary tremendously. For instance, for 
congelation ice frozen under calm conditions, the crystals are 
vertically elongated columns, and each crystal itself is a 
composite of pure ice platelets separating layers of brine 
inclusions. The orientation of each individual crystallite 
is determined by the direction of the 
c-axis, which is perpendicular 
to the orientation of platelets or lamellae of pure ice.
Finding the bounds for the crystal orientations enables us to electromagnetically distinguish columnar ice from granular ice. This is a critical problem in sea ice physics and biology, as these different structures have vastly different fluid flow properties, % (with 5\% \textit{vs.} 10\% brine volume fraction at the percolation threshold), 
which affects melt pond evolution, nutrient replenishment, brine convection, and other mesoscale processes in the ice cover. 

%an algorithm of estimating the effective angle of the deviation of the polycrystal's axis from the vertical direction. 

Bounds for the effective permittivity of polycrystalline composites are much tighter than those bounding the permittivity of a general two-component material. Such polycrystalline bounds constructed in \citep{Gully:PRSA:2015}
are shown in Fig. \ref{inv-bounds}(c). This dramatic improvement over the classic two-component bounds is due to additional information about single crystal orientations included in the new bounds.  %Polycrystalline bounds  for the permittivity sea ice (left)  \citep{Gully:PRSA:2015} 
%(with the measured data on permittivity of sea ice \citep{Arcone:JGR-14281})
%provide a much tighter bound than general two-component material and statistically isotropic two-component material for sea ice  given on the right (notice a different scale). This dramatic improvement over the classic two-component bounds is due to additional information about single crystal orientations included in the new bounds. 

As was discussed in the polycrystal section, 
the zeroth moment $ \mu_{kk}^0$ in  \eqref{eq:Mass_mu} of the measure $\mu$ in the integral representation of the effective properties of a uniaxial polycrystalline material  is 
%
$ \mu_{kk}^0=\langle |X_1\vece_k|^2\rangle $, 
%
where the statistical average $\langle |X_1\vece_k|^2\rangle$ %in~\eqref{eq:Mass_mu}
can be viewed as the ``mean crystal orientation''  related to the percentage of
crystallites oriented in the $k^{\text{th}}$ direction. Extending the inverse bounds method \citep{Cherkaev:1996,Cherkaev:WRM-437,Tripp:1998}
to polycrystalline materials, the inverse polycrystalline bounds \citep{Gully:PRSA:2015} estimate the mean crystal orientation by bounding the zeroth moment  $ \mu_{kk}^0$ of the measure $\mu$ using measured data on the ice  permittivity. This procedure gives an analytic estimate (the first order inverse bounds) for the range of values of the mean crystal orientation similar to (\ref{inv1b}): 
\begin{eqnarray}                                 \label{eq:invpoly}
\langle e_k^T X_1  e_k \rangle_l \leq \langle e_k^T X_1  e_k \rangle \leq \langle e_k^T X_1  e_k \rangle_u\,, \qquad \qquad \qquad     \nonumber \\ 
\langle e_k^T X_1  e_k \rangle_l = |f|^2 \frac{Im(\overline{s})}{Im(f)}\,,
\qquad                              
\langle e_k^T X_1  e_k \rangle_u = 1 -  |g|^2   \frac{Im(\overline{t}) }{Im(g)}\,,
\end{eqnarray}
Here, $X_1$ is defined in the polycrystalline section as $X_1=R^{\,T}CR$, 
$f$ is the known value of $F(s)$ 
and $g$ is the known value of $G(t)=1-\epsilon^*/\epsilon_1$  with $t=1-s$.

Inverse polycrystalline bounds computed in \citep{Gully:PRSA:2015} for different types of sea ice, granular and columnar, show that the method allows the type of ice to be revealed through electromagnetic data. 
%
For statistically isotropic granular ice shown in Fig. \ref{polycrystalline}(a)-top, the inverse mean crystal orientation bounds \citep{Gully:PRSA:2015} estimate 
the deviation angle as $\pi/2 \pm .02$ (with the true value $\pi/2$).
%
The inverse mean crystal orientation bounds \citep{Gully:PRSA:2015} for columnar ice 
(see Fig. \ref{polycrystalline}(a)-bottom), 
estimate  the angle of deviation of the crystal's axis from the vertical as $20^o \pm 8^o.$ 
These results demonstrate a significant difference in the reconstructed mean orientations of crystals in columnar and in granular ice and provide a foundation for distinguishing the types of ice using electromagnetic measurements. 


%*******  Generalizations: Kenzie   
%Generalization of these polycrystalline bounds to the case when the $c$-axis has a  
%Gaussian distribution with known mean angle and variance in the horizontal plane is developed in \citep{McLean:poly:2019} as a method for obtaining bounds on effective permittivity of columnar sea ice that has a preferred direction in the horizontal plane due to a prevailing ocean current.

% END INVERSE 







\section{Analytic continuation for advection diffusion processes}
The enhancement of diffusive transport of passive scalars by complex
fluid flow plays a key role in many important processes in the global 
climate system~\citep{Washington:1986:9780935702521} and Earth's
ecosystems~\citep{Lorenzo:2013:26:4}. Advection of geophysical
fluids intensifies the dispersion and large scale transport of
heat~\citep{Moffatt:RPP:621},
pollutants~\citep{Bilger:10.1175,Beychok:1994:9780964458802,Samson:1988:88009978},
and
nutrients~\citep{Lorenzo:2013:26:4,Hofmann:ANS:2004:265075} diffusing
in their environment. 
% Advective processes also enhance the large scale
% transport of plankton~\citep{Hofmann:ANS:2004:265075}, which is an
% important component of the food web that sustains life in the polar
% oceans.  
% The transport of vast sea ice floes in the polar oceans is
% driven by a seasonally and regionally changing balance in oceanic and
% atmospheric
% forces~\citep{Lukovich:2011:JGRC12087,Washington:1986:9780935702521}. These
% forces
% %are vorticity dominated and
% can enhance the transport of sea
% ice floes by eddie fluxes in the atmosphere and ocean
% currents~\citep{Watanabe:2009JPO4010,Rampal:2008:JGRC10653,Rampal:2009:JGRC11312,Lukovich:2011:JGRC12087}. 
% %as well as intense jets through straits and
% %canyons~\citep{Watanabe:2009JPO4010}.  
%
In sea ice dynamics,
where the ice cover couples the atmosphere to the polar
oceans \citep{Washington:1986:9780935702521}, the transport of sea 
ice can also be enhanced by eddy
fluxes and large scale coherent structures in the ocean
\citep{Watanabe:2009JPO4010,Lukovich:AG:2015,Dinh:anamolous:2021}. 
In sea ice thermodynamics, the temperature field of the 
atmosphere is coupled to the temperature field of the ocean through
sea ice, a composite of pure ice with brine inclusions whose
volume fraction and connectedness depend strongly on temperature
\citep{Thomas:2003:SI,Golden:GRL:L16501,Golden:NAMS:2009}. Convective brine
flow through the porous microstructure can enhance thermal transport
through the sea ice layer
\citep{Lytle:JGR-8853,Worster:PTRSA:2015} %,Kraitzman_thermal_cond_2023}. %,Liu:2015}. 

% Above, we discussed the analytic continuation method (ACM) for representing 
% transport in composite materials. This method reduces the analysis of 
% complex composite materials, with rapidly varying structures in space, 
% to solving averaged, or \emph{homogenized} equations that do not have 
% rapidly varying data, and involve an effective parameter. We also 
% discussed how the ACM leads to a \emph{Stieltjes integral 
% representation} of the effective parameters for composites, involving
% the spectral measure of a self-adjoint operator (or matrix). This 
% integral representation yields rigorous \emph{forward bounds} for the 
% effective parameters of composites, given partial information on the 
% microgeometry \citep{Bergman:PRL-1285,Milton:APL-300,Golden:CMP-473,Bergman:AP-78}.
% One can also use the integral representations to obtain \emph{inverse
%   bounds}, allowing one to use data about the electromagnetic response
% of a sample, for example, to bound its structural parameters, such as
% the volume fraction of each of the components \citep{Golden:JBM:337}.  


% The effective parameter problem for two-component composites discussed in Section \ref{sec:two_component} reduces the 
% analysis of complex composite materials, with rapidly varying structures in 
% space, to solving averaged, or \emph{homogenized} equations that do not have 
% rapidly varying data, but involve an effective parameter.
% This mathematical framework was extended 
% \citep{McLaughlin:SIAM_JAM:780,Majda:Kramer:1999:book}
% to homogenize the advection diffusion equation, involving a complex fluid
% velocity field with rapidly varying structures in both space and time, yielding 
% a homogenized equation involving an effective diffusivity tensor $\Dm^*$.
% Moreover, a Stieltjes integral representation for $\Dm^*$ was obtained
% \citep{Avellaneda:PRL-753,Avellaneda:CMP-339,Murphy:ADSRF-2019,Murphy:ADSTPF-2017,Pavliotis:PHD_Thesis,Avellaneda:PRE:3249,Bhattacharya:AAP:1999:951}, 
% which involves a spectral measure $\nu$ of a self-adjoint, Hermitian operator 
% (or matrix). This representation separates the molecular diffusivity 
% %$\varepsilon$ 
% %(or equivalently the P\'{e}clet number $Pe$) 
% (or equivalently the P\'{e}clet number) 
% from the geometric complexity and dynamics 
% of the flow \citep{Murphy:ADSRF-2019,Murphy:ADSTPF-2017}, 
% which are incorporated through 
% the moments $\nu_n=\int\lambda^n \,d\nu(\lambda)$ of the measure $\nu$. 
% As in the setting of composite media, Pad{\'e} approximants provide
% rigorous bounds on the components of $\Dm^*$
% involving the moments $\nu_n$ \citep{Baker:1996:Book:Pade,Avellaneda:CMP-339}. 
% The bounds become progressively tighter as additional moments are 
% incorporated. 







% Over the years a broad range of mathematical techniques have been
% developed that reduce the analysis of complex composite materials, with
% rapidly varying structures in space, to solving averaged, or
% \emph{homogenized} equations that do not have rapidly varying data,
% and involve an effective parameter. The effective parameter problem
% regarding multicomponent composite media for thermal conductivity, or
% equivalently \citep{MILTON:2002:TC} electrical conductivity and
% permittivity, and magnetic permeability, was developed
% \citep{Papanicolaou:RF-835} by Papanicolaou and Varadhan for the
% (steady state) static limit. This multicomponent work was
%  extended \citep{Golden:JSP-655} by Golden and Papanicolaou to the
% frequency dependent quasi--static case (with complex parameters).  
% %The full time dependent problem for the effective dielectric
% %coefficient is given in \citep{Simeonova:MMS:1113}.

% The effective parameter problem regarding \emph{two--component} composite media, in
% the quasi--static limit, was developed \citep{Golden:CMP-473} by Golden and Papanicolaou, leading to a \emph{Stieltjes integral
% representation} of the effective parameters for composites,
% discussed in Section \ref{sec:two_component}. This
% analytic continuation method (ACM) is based 
% on the spectral theorem and a resolvent formula involving a random 
% self--adjoint operator (or matrix). This integral representation
% yields rigorous \emph{forward bounds} for the effective parameters of
% composites, given partial information on the microgeometry
% \citep{Bergman:PRL-1285,Milton:APL-300,Golden:CMP-473,Bergman:AP-78}.
% One can also use the integral representations to obtain \emph{inverse
%   bounds}, allowing one to use data about the electromagnetic response
% of a sample, for example, to bound its structural parameters, such as
% the volume fraction of each of the components \citep{Golden:JBM:337}.  



%% copied from above
Over the years, a broad range of mathematical techniques have been
developed that reduce the analysis of complex composite materials with
rapidly varying structures in space, to solving averaged or
\emph{homogenized} equations that do not have rapidly varying coefficients and involve an effective parameter.
%
The core idea here is that the motion of a particle diffusing 
in a velocity field with regular geometry (including stationary random, periodic, or quasiperiodic variations)
is governed on a large scale and long times by a diffusion equation  % akin to Brownian motion
with an effective diffusion tensor which depends on the geometry of the 
velocity field and the local diffusivity.
In \citep{Taylor:PRSL:196}, it was first shown that a long time, 
large scale dispersion of passive scalars can be described by an effective diffusivity
tensor $\Dm^*$.
Motivated by \citep{Papanicolaou:RF-835}, the effective parameter
problem was extended to complex velocity fields, with rapidly varying
structures in both space and time, providing a rigorous mathematical foundation 
for calculating effective (eddy)
viscosity and the effective (eddy) diffusivity tensors
\citep{McLaughlin:SIAM_JAM:780}. 
The effective parameter problem of
(anomalous) super--diffusion and sub--diffusion  is given in
\citep{Biferale:PF:2725,Fannjiang:PD:145}. Based on
\citep{McLaughlin:SIAM_JAM:780}, %Avellaneda and Majda
\cite{Avellaneda:PRL-753,Avellaneda:CMP-339} adapted the ACM
\citep{Golden:CMP-473} to the advection diffusion equation and obtained
a \emph{Stieltjes integral representation} of the effective diffusivity
tensor $\Dm^*$ for flows with zero mean drift involving the
P\'{e}clet number $\xi$ of the flow. This representation
encapsulates the geometric complexity of the flow in a spectral
measure associated with a random Hermitian operator (or
matrix). Mimicking methods developed for composite media 
\citep{MILTON:2002:TC}, they obtained rigorous bounds on the components
of $\Dm^*$. Moreover, in analogy to methods developed for
composites \citep{MILTON:2002:TC}, they also found velocity fields which
realize these bounds, such as the famous confocal sphere
configurations that realize the Hashin--Shtrikman bounds of composites
%media
\citep{Hashin:JAP-3125,Avellaneda:CMP-339}. Remarkably, this 
method has also been extended to time dependent flows
\citep{Avellaneda:PRE:3249,Murphy:ADSTPF-2017}, flows with incompressible \emph{nonzero}
effective drift \citep{Pavliotis:PHD_Thesis,Fannjiang:1994:SIAM_JAM:333},
flows where particles diffuse according to linear collisions
\citep{Pavliotis:IMAJAM:951}, and solute transport in porous media
\citep{Bhattacharya:AAP:1999:951}, which has a direct 
application to diffusive brine advection in sea ice. All yield Stieltjes
integral representations of the symmetric and, when appropriate, the 
antisymmetric part of $\Dm^*$.



We now briefly describe our recent %mathematical 
results on
this framework 
% in terms of the effective 
% thermal conductivity of sea ice in the presence of a brine velocity field
%\citep{Kraitzman_thermal_cond_2023,
\citep{Murphy:ADSTPF-2017,Murphy:ADSRF-2019}.
It is an % well defined and 
important example of how Stieltjes integral
representations can provide a rigorous basis for analysis
of problems for sea ice involving advection diffusion processes.
% This framework can also be used to model the advection 
% enhanced dispersion of sea ices floes by winds 
% and ocean currents collectively giving rise 
% to a transporting flow, which we will discuss below
% in our proposed work.
% While this proposal focuses primarily on broadening the
% applicability of homogenization techniques to larger scales,
% we also have a few ongoing smaller scale homogenization
% projects such as moment calculations in the advection
% diffusion framework, and work on the thermal conductivity
% of sea ice, the electromagnetic properties of polycrystalline
% media, and the evolution of melt ponds on Arctic sea ice.
% 
% 
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sea ice thermal conductivity %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The enhancement of sea ice thermal conductivity %$\kappa$ 
% above its value $k_0$ with no fluid flow,
% by a space-time periodic brine velocity field $u$,
% is described by the advection diffusion equation 
% \citep{Lytle:JGR-8853,Worster:PTRSA:2015,Kraitzman_thermal_cond_2023,Murphy:ADSRF-2019,Murphy:ADSTPF-2017}
% %
% \begin{align}\label{eq:ADE}
%   \frac{\partial T}{\partial t} 
%   +\phi \,u \cdot\nabla T\,
%   =\frac{k_0}{\rho\,c}\,\Delta T,
%   \qquad
%   T(0,x)=T_0(x).
%   % \quad
% %   t>0,
% %   \quad
% %   \vecx\in\mathbb{R}^d.
% % \qquad
% % \bnabla\bcdot\vecu=0.
% \end{align}
% %
% %for $t>0$ and $\vecx\in\mathbb{R}^d$.
% Here, we take $k_0=2$ W/(m W), $\rho=960\,$kg/m$^3$
% is the bulk density of sea ice, 
% $c= 2.11$ J/(g K) is the specific heat,  and $\phi$
% is the brine volume fraction, a function of temperature and salinity, 
% though taken here as an independent 
% parameter for simplicity. Moreover, the 
% incompressible
% \emph{bulk} brine velocity field 
% $u(t,x)$ and the initial temperature $T_0(x)$ are given. 
% %assumed 
% %given, $u$ is incompressible, 
% %$\nabla\cdot u=0$. 
% %and $\partial_t$ denotes differentiation with respect to time $t$.
%
%%%%%%%%%%%%%%%%%%%%%%%%%
% Effective diffusivity %
%%%%%%%%%%%%%%%%%%%%%%%%%
% The enhanced diffusive transport of temperature $T(t,\vecx)$ through sea brine 
% microstructure, above the bare molecular diffusivity value $\kappa_0$, 
% by a space-time periodic brine velocity field $\vecu(t,\vecx)$ 
% can be described by the advection-diffusion equation       
% 
The dispersion of a cloud of passive scalars with density $\phi(t,x)$
(not to be confused with earlier usage as the brine volume fraction)
diffusing with molecular diffusivity $\varepsilon$ and being advected by a
%random,
incompressible velocity field $\vecu(t,x)$ satisfying
$\bnabla\bcdot\vecu=0$ is described by the advection-diffusion
equation         
\begin{align}\label{eq:ADE}
\frac{\partial \phi}{\partial t} 
= \,u \cdot\nabla \phi\,
+\epsilon\,\Delta \phi,
\qquad
\phi(0,x)=\phi_0(x).
% \quad
%   t>0,
%   \quad
%   \vecx\in\mathbb{R}^d.
% \qquad
% \bnabla\bcdot\vecu=0.
\end{align}
%
%for $t>0$ and $\vecx\in\mathbb{R}^d$.
Here, the initial density $\phi_0(\vecx)$ and the fluid velocity field
$\vecu$ are assumed to be given.
% , and $\vecu$ satisfies $\bnabla\bcdot\vecu=0$. 
In equation~\eqref{eq:ADE}, the molecular diffusion constant is $\varepsilon>0$,
% $d$ is the spatial dimension of the 
% system, 
% $\partial_t$ denotes partial differentiation with respect to time
% $t$, 
and $\Delta=\bnabla\bcdot\bnabla =\nabla^2$ is the Laplacian. 
%
This equation also models the transport of heat advected by the fluid velocity field $u$ and diffused with molecular diffusion 
coefficient $\varepsilon$.
%
%
%Moreover,$\vecpsi\bcdot\vecvarphi=\vecpsi^{\,T}\overline{\vecvarphi}$, $\vecpsi^{\,T}$ denotes transposition of the vector $\vecpsi$, and $\overline{\vecvarphi}$ denotes component-wise complex conjugation, with $\vecpsi\bcdot\vecpsi=|\vecpsi|^2$. Later, we will use this form of the dot product over complex fields, with built in complex conjugation. However, we emphasize that all quantities considered in this section are \emph{real-valued}. The random paths of the tracer particles are determined~\citep{Fannjiang:1997:1033} by the stochastic differential equation   
%
%\begin{align}
%\d\vecx(t)=\vecu(t,\vecx(t))\d t+\sqrt{2\varepsilon}\;\d W(t),
%\quad
%\vecx(0)=\vecx_0,
%\end{align}
%
%with the initial tracked tracer particle location $\vecx_0$ given and 
%$W(t)$ is standard Brownian motion (the Wiener process).
%
% Moreover, it is assumed~\citep{McLaughlin:SIAM_JAM:780} that
% the fluid velocity field $\vecv(\vecx)$ is stationary and that the
% correlation functions $R_{ij}(z)=\langle v_i(\vex+\vecz)\,v_j(\vecx)\rangle$ have a
% Fourier transform that is continuous at the origin.
%
%
To simplify our presentation, we assume that the velocity field $u$ in equation \eqref{eq:ADE} is 
temporally and spatially periodic. 
Non-dimensionalization \citep{Murphy:ADSRF-2019} and homogenization \citep{McLaughlin:SIAM_JAM:780} of  
equation \eqref{eq:ADE} shows that long time, 
large volume (or area) \emph{macroscopic} thermal transport is 
described by a diffusion equation involving an 
averaged scalar density $\bar{\phi}$ and
% averaged temperature $\bar{T}$ and
a symmetric, constant \citep{Pavliotis:PHD_Thesis}
effective diffusivity % conductivity 
tensor $\bkappa^*$ \citep{McLaughlin:SIAM_JAM:780},  
%
% The long time, large scale dispersion of the passive scalars can be
% described~\citep{Taylor:PRSL:196} by an effective diffusivity
% tensor $\Dm^*$. An explicit representation for $\Dm^*$ can be found
% using methods of homogenization
% theory~\citep{Bensoussan:Book:1978,Papanicolaou:RF-835}. These methods
% have demonstrated that the averaged or \emph{homogenized} behavior of
% the advection-diffusion equation in~\eqref{eq:ADE} is determined by a
% diffusion equation involving an averaged scalar density $\bar{\phi}$ and
% a (constant) effective diffusivity tensor 
% $\Dm^*$~\citep{McLaughlin:SIAM_JAM:780}   
%
\begin{align}\label{eq:phi_bar}
\frac{\partial \bar{\phi}(t,x)}{\partial t}
%\partial_t\,\bar{T}(t,x)
=\nabla\cdot[\bkappa^*\,\nabla \bar{\phi}(t,x)], \quad
\bar{\phi}(0,x)=\phi_0(x).
\end{align}
%
% end blue



For simplicity, we consider a diagonal coefficient $\kappa_{kk}^*$, $k=1,\ldots,d$, 
of $\bkappa^*$, set $\kappa^*=(\bkappa^*)_{kk}$, 
and write $u=u_0v$, where $u_0$ is the (constant) strength of $u$ and $v$ is a 
non-dimensional velocity field containing the geometric and dynamic information about $u$.
% and write the non-dimensional velocity field
% $v=(\phi\,\ell/\tau)\,u$ and sea ice thermal ``diffusivity''
% $\varepsilon=k_0\tau/(\ell^2\,\rho\,c)$, where $\tau$ and $\ell$ 
% are system time and length scales, then 
% \citep{McLaughlin:SIAM_JAM:780,Avellaneda:CMP-339,Murphy:ADSRF-2019,Murphy:ADSTPF-2017} 
% The components $\Dmc^*_{jk}$, $j,k=1,\ldots,d$, of the 
% tensor $\Dm^*$ are given
% by~\citep{McLaughlin:SIAM_JAM:780}   
%
%
In these non-dimensional variables the P\'eclet number $\xi$ and molecular diffusivity
$\varepsilon$ are related by $\xi=1/\varepsilon$ \citep{Murphy:ADSTPF-2017,Murphy:ADSRF-2019}.

Using a mathematical framework that is strikingly similar to that in 
Section \ref{sec:two_component}, the effective diffusivity has the following
Stieltjes integral representation 
\citep{McLaughlin:SIAM_JAM:780,Avellaneda:CMP-339,Murphy:ADSTPF-2017,Murphy:ADSRF-2019}
% involving the spectral measure $\nu$ of a random operator $i\Gamma S\Gamma$
%
\begin{align}\label{eq:Eff_Diffusivity_Sobolev}
%D^*_{jk}=\varepsilon(\delta_{jk}+\langle w_j,\chi_k\rangle_{1,2}),
\kappa^*=\varepsilon(1+\langle|\nabla w_k|^2\rangle),
\qquad
\langle|\nabla w_k|^2\rangle=\int_{-\infty}^\infty\frac{\d\nu(\lambda)}
{\varepsilon^2+\lambda^2}\,,
%   % w_k=(\varepsilon+A)^{-1}g_k, \qquad 
%   \chi_k=(\varepsilon+A)^{-1}(-\Delta)^{-1}(\beta \,v_k), 
%   \qquad 
%   %A=\Delta^{-1}(\vecu \bcdot\bnabla -\partial_t), \quad
%   %g_k=(-\Delta)^{-1}u_k.
%   %\Am^*_{jk}=\langle A w_j,w_k\rangle_{1,2}\,,
%   %\quad
%   A=(-\Delta)^{-1}(\partial_t+\beta \,v \cdot\nabla)\,.
\end{align}
%
where $\langle\cdot\rangle$ denotes averaging over the space-time period cell for periodic flows \citep{Murphy:ADSTPF-2017,Murphy:ADSRF-2019} or statistical
average for random flows \citep{Avellaneda:PRL-753,Avellaneda:PRE:3249}
and $w_k$ is the solution to a \emph{cell problem} 
\citep{McLaughlin:SIAM_JAM:780,Murphy:ADSTPF-2017}.
An equivalent statement which emphasizes the connection to the two component 
composites setting in equation \eqref{eq:F(s)} is
%
\begin{equation}
F(\varepsilon) = 1 -  \frac{\kappa^*}{\varepsilon} = 
\int_{-\infty}^\infty\frac{\d\nu(\lambda)}
{\varepsilon^2+\lambda^2}\,.
\label{F_for_diffusivity}
\end{equation}
%
The vector field $E(t,x)=\nabla w_k(t,x)+e_k$ satisfies equation
\eqref{eq:Maxwells_Equations} for two-component composite materials, with $D=\beps E$, 
$\beps=\varepsilon I + S$, $S=(-\Delta)^{-1}\partial_t+H$, and $\beps$ plays the role
of the medium's electrical permittivity tensor \citep{Murphy:ADSTPF-2017,Murphy:ADSRF-2019}. 
Here, $\partial_t$ denotes partial differentiaion with respect to time and
$H(t,x)$ is the \emph{stream matrix}, 
given in terms of the incompressible velocity field $v=\nabla\cdot H$ and satisfies
$H^T=-H$ \citep{Avellaneda:CMP-339,Avellaneda:PRL-753}. When the flow is time-independent,
$v=v(x)$, then $ w_k= w_k(x)$ and $S=H(x)$. Moreover 
$\kappa^*=\epsilon^*$, with $\epsilon^*=(\beps^*)_{kk}$ 
defined above \citep{Murphy:ADSTPF-2017}. % in equation \eqref{eq:sigma*}. 
The integral representation for $\kappa^*$ in Equation
\eqref{eq:Eff_Diffusivity_Sobolev} follows from the resolvent formula
%$\nabla w_k=(\varepsilon I + \imath M)[-\Gamma H e_k]$, 
% $\nabla w_k=(\varepsilon I + i \Gamma S\Gamma)^{-1}[-\Gamma H e_k]$,
%
\begin{align}\label{eq:resolvent_ADE}
\nabla w_k=(\varepsilon I + \Gamma S\Gamma)^{-1}g_k,
\quad
g_k=-\Gamma H e_k,
\end{align}
%
which is an analogue of Equation
\eqref{eq:Electric_Field_two_component}. 
The operator $\Gamma S\Gamma$ is \emph{antisymmetric}
due to the asymmetry of both the operators $\partial_t$ and $H$, so $i
\Gamma S\Gamma$ is a \emph{self-adjoint} operator \citep{Murphy:ADSTPF-2017}, 
where $i=\sqrt{-1}$ is the imaginary unit
and $\Gamma=-\nabla(-\Delta)^{-1}\nabla\cdot\,$ is the same projection operator 
arising in the setting of two-component composites.
Equation \eqref{eq:Eff_Diffusivity_Sobolev} shows that 
%the effective diffusivity in the principle directions $\Dmc^*_{kk}$, are 
brine advection \emph{enhances} the thermal 
diffusivity 
%(and the thermal conductivity) of sea ice,
%$D^*$ is \emph{enhanced} above the bare molecular diffusivity value 
%$\Dmc^*=\varepsilon(1+\langle|\bnabla\chi_k|^2\rangle)\ge\varepsilon$.  
since $\kappa^*\ge\varepsilon$. 
%of sea ice, as thermal diffusivity and conductivity are directly proportional to one another. 

% Here, 
% %$\langle f,h\rangle_{1,2}=\langle\bnabla f\bcdot\bnabla h\rangle$ is a Sobolev-type
% %\emph{sesquilinear} 
% %inner-product \citep{Murphy:ADSTPF-2016},
% $\langle\cdot\rangle$ denotes averaging over the space-time period cell, $v_k$ is the 
% $k$th component of $v$,
% and $A$ is a skew-adjoint operator which can be written as $A=\imath M$,
% where $M$ is \emph{self-adjoint} \citep{Murphy:ADSTPF-2016}. 
% Equation \eqref{eq:Eff_Diffusivity_Sobolev} shows that 
% %the effective diffusivity in the principle directions $\Dmc^*_{kk}$, are 
% brine advection \emph{enhances} the thermal conductivity of sea ice
% %$D^*$ is \emph{enhanced} above the bare molecular diffusivity value 
% %$\Dmc^*=\varepsilon(1+\langle|\bnabla\chi_k|^2\rangle)\ge\varepsilon$.  
% $\kappa^*\ge\varepsilon$. The spectral theorem then yields a Stieltjes 
% integral representation for $\kappa^*$ \citep{Murphy:ADSTPF-2016}
% %
% \begin{align}\label{eq:Integral_Rep_kappa*}
% %   D^*_{jk}=\varepsilon\left(\delta_{jk}+\int_{-\infty}^\infty\frac{\d\nu_{jk}(\lambda)}{\varepsilon^2+\lambda^2}\right),
%   \kappa^*=\varepsilon(1+S(\varepsilon)),
%   \quad
%   S(\varepsilon)=\int_{-\infty}^\infty\frac{\d\nu(\lambda)}
%                                                 {\varepsilon^2+\lambda^2}\,,
% %\Am^*_{jk}=-\int_{-\infty}^\infty\frac{\lambda\,\d\Imag\,\mu_{jk}(\lambda)}{\varepsilon^2+\lambda^2}\,.    
% \end{align}
% %
% where $\nu$ is a positive \emph{spectral measure} associated with the 
% self-adjoint operator $M=-\imath A$.


% Writing the Stieltjes integral $S(\varepsilon)$ as a \emph{series of Stieltjes} in the 
% variable $z=\varepsilon^{-2}$ 




\begin{figure}%[h!] 
\centering
\includegraphics[width=\textwidth]{BCFlow_spectral_measures.png}
\caption{\textbf{Spectral measures for effective diffusivities.}
(a) Streamlines for BC-flow with fluid velocity field $v=(C\cos{y},B\cos{x})$ with $B=1$, $C=1-p$, and increasing values of $p$ from left to right. BC-cell-flow with closed streamlines arises for $B=C=1$. As the value of $C$ decreases from 1, the streamlines elongate in the y-direction giving rise to large scale thermal transport even in the absence of molecular diffusivity, when $\varepsilon=0$. BC-shear flow is attained when $B=1$ and $C=0$, for which the spectral measure is known to be a $\delta$-function at the origin $\lambda=0$ \citep{Avellaneda:CMP-339}, as shown in the rightmost panel of (a). (b) Corresponding spectral functions $\nu_{22}(\lambda)$ (histogram representations) for the spectral measures $\nu_{22}$ in (c), which display the spectral measure weights $m_i^2$ vs eigenvalues $\lambda_i$ of the matrix $G=i\Gamma H \Gamma$. Since the streamlines do not become elongated in the x-direction, the spectral functions for the measure $\nu_{11}$ are qualitatively similar to that of $\nu_{22}$ for BC-shear flow, for all $p$. The spectral functions in (b) are ensemble averaged for $B=1$ and $C=1-\zeta$ with $\zeta\sim U([0,p])$, with $p=0.01$, 0.5, 0.7, and 1, from left to right. The zoomed in insets in (c) with $-5 \times 10^{-4} \le \lambda_i \le 5\times 10^{-4}$ show the density of the measure near the spectral origin $\lambda=0$ increasing as BC-cell flow transitions to BC-shear flow.  } 






%(a) Streamlines for BC-flow with velocity field $v=(C\cos{y},B\cos{x})$ and $B=C=1$. (b) Pad{\'e} approximant upper $\kappa^*[M/M]$ and lower $\kappa^*[M-1/M]$  bounds for $\kappa^*$, for various values of $M$, calculated %via equation \eqref{eq:Pade_ADE} 
%for BC-flow with $C=B$, as a function of the flow strength $B$. (c) The spectral function (spectral masses $m_j$ versus eigenvalues $\lambda_j$) computed via analogues of equations \eqref{eq:Discrete_E_F(s)} and \eqref{eq:Stieltjes_F_Discrete_polycrystal} \citep{Murphy:ADSTPF-2017}.

\label{fig:BC-Flow}      
\end{figure}




Analytical calculations of the spectral measure $\nu$ are extremely difficult except for simple flows like
shear flow \citep{Avellaneda:CMP-339}. However, 
Pad{\'e} approximants $[L/M]$
provide rigorous, converging 
upper and lower
bounds \citep{Baker:1996:Book:Pade}
for the \emph{Stieltjes function}
$f(z)=\langle|\nabla w_k|^2\rangle/z=F/z$ in 
equations
\eqref{eq:Eff_Diffusivity_Sobolev} and \eqref{F_for_diffusivity}; 
Pad{\'e} approximants can be calculated using the moments $\nu_n$ of $\nu$. 
%in the variable 
%With $z=\varepsilon^{-2}$ and using the moments $\nu_n$
%$\nu_n=\int\lambda\,\d\mu(\lambda)$ 
% of the measure  $\nu$, % one achieves the characterization
%which is associated with the P{\'e}clet number $Pe$
% by writing it as a \emph{series of Stieltjes}
% $S(z)=z\sum_{n=0}^\infty(-1)^n\,\nu_n\,z^n$ 
%
%\begin{align}\label{eq:Pade_ADE}
%$[M-1/M]\le f(z) \le [M/M]$,
%\qquad
%$f(z)=\sum_{n=0}^\infty(-1)^n\,\nu_{\,2n}\,z^n$.
% \quad
% \,,
% \qquad
% \nu_n = \int_{-\infty}^\infty\lambda^n \d\nu(\lambda).
%\end{align}
%
%where %$M-1,M=0,1,2,\ldots$ and 
%the summation follows from writing $1/(1+z\lambda^2)$ as a 
%geometric series. Bounds for $\kappa^*$ can also be obtained using variational methods 
%\citep{Avellaneda:CMP-339,Fannjiang:1994:SIAM_JAM:333}. While the mathematical framework 
%summarized by equations \eqref{eq:phi_bar}--\eqref{eq:Pade_ADE} provides an elegant,
%analytic representation for $\kappa^*$ and a means to obtain tight bounds on $\kappa^*$, 
Despite this, the lack of a method to calculate 
the moments $\nu_n$ of $\nu$ has impeded progress on 
obtaining explicit 
bounds for specific flows using this %bounding 
procedure 
\citep{Avellaneda:CMP-339,Avellaneda:PRL-753} since 1991! 
We have recently developed a mathematical framework %\citep{Murphy_Pade_Bounds} 
that can be used to compute, in principle, \emph{all} of the moments $\nu_n$ for a
spatially or space-time periodic velocity field $v$, hence Pad{\'e} approximant bounds. 
We have utilized these results to provide rigorous bounds for the enhancement of sea ice thermal conductivity by brine fluid velocity fields. These findings will be published elsewhere. 

%\citep{Kraitzman_thermal_cond_2023}.
%
% Results for BC-flow, with $v=(C\,\cos{y},B\,\cos{x})$ and $B=C$ are 
% shown in \figref{fig:BC-Flow}(c). 
% With very little existing theoretical 
% work on the effective thermal conductivity of sea ice with a non-zero
% brine velocity field, refining our bounds
% to capture various field data sets is ongoing
% \citep{Kraitzman_thermal_cond_2023,Murphy:ADSRF-2019}.


%This framework is equivalent to 
%that described above, reproducing the result in 
%equation
%\eqref{eq:Eff_Diffusivity_Sobolev} 
%with the same measure $\nu$ \citep{Murphy:ADSTPF-2017,Murphy:ADSRF-2019}. Instead of focusing 
%on the vector field $\nabla w_k$, the scalar field $ w_k$ is used, which has the resolvent 
%formula $ w_k = (\epsilon + G)^{-1}[(-\Delta)^{-1}v_k]$, $G=\imath (-\Delta)^{-1}D_t$ is a 
%self-adjoint operator, 
%on an appropriate Hilbert space, 
%$v_k$ is the $k$th component of $v$, $D_t=\partial_t+v\cdot\nabla$ 
%for time-dependent $v$ and $D_t=v\cdot\nabla$ for time-independent $v$. Using this approach, 
%the moments $\nu_n$ of $\nu$, hence the bounds in equation \eqref{eq:Pade_ADE} can be 
%calculated exactly \emph{in closed form}. Our recent results for BC-flow, with 
%$v=(C\,\cos{y},B\,\cos{x})$ and $B=C$ are displayed in Figure 
%\ref{fig:BC-Flow}.  

%We propose to further develop and investigate rigorous
%bounds on the convection enhanced thermal conductivity 
%of sea ice with brine velocity fields. Field measurements
%of this key parameter are notoriously difficult, and
%there is very little theoretical work in this direction.
%We will also further develop this framework and apply it to a 
%broad range of 2D and 3D steady periodic flows, 
%as well as 2D and 3D space-time periodic flows, 
%as well as spatially periodic, time-stochastic flows. 
%We also propose to develop inverse 
%bounds
%\citep{McPhedran:APA:19,McPhedran:MRSP:1990:195,Cherkaev:WRM-437,Cherkaev:IP-1203,Cherkaev:IP-065008,Zhang:JCP-5390,Bonifasi-Lista:PMB-3063,Cherkaev:JBiomech-345,Day:1999,Golden:JBM:337}, 
%which incorporate measurements of the effective thermal 
%conductivity of sea ice 
%the effective diffusivity of sea ice floes, via satelite bouy data, for example, 
%to potentially recover geometric and dynamic information 
%about the underlying fluid velocity field. 


%\subsection{Spectral measure computations for advection diffusion processes} 
%\label{ADE_discrete}
%


\vspace{3ex}


\noindent
\textbf{Spectral measure computations for advection diffusion processes.} 
We have extended our numerical methods 
discussed for two component media in Section \ref{sec:two_component_discrete} 
to compute the spectral measure $\nu$ for 
spatially periodic flows \citep{Murphy:ADSRF-2019} and developed Fourier methods
for computing $\nu$ for space-time periodic flows \citep{Murphy:ADSTPF-2017}. 
For simplicity, we focus our discussion to the setting of spatially periodic flows.
Computing the spectral
measure $\nu$ for a given flow involves discretizing the spatially dependent
stream matrix $H(x)$, which becomes a banded antisymmetric matrix satisfying
$H^T=-H$. The projection matrix $\Gamma$ is given by that in Section
\ref{sec:two_component_discrete} and the key self-adjoint operator is given
by $G=i\Gamma H\Gamma$, which becomes a \emph{Hermitian} matrix. 
In this
case, the integral in \eqref{eq:Eff_Diffusivity_Sobolev} and the resolvent
in \eqref{eq:resolvent_ADE} are given in terms of the eigenvalues $\lambda_i$ and 
eigenvectors $\vecw_i$ of the matrix $G$
%
\begin{align}
\nabla w_k
=\sum_i\frac{m_i}{\varepsilon+i\lambda_i}\vecw_i\,,
\quad
\langle|\nabla w_k|^2\rangle
=\sum_i\left\langle\frac{|m_i|^2}{\varepsilon^2+\lambda_i^2}\right\rangle\,,
\quad
m_i=\vecw_i\cdot g_i\,,
\end{align}
%
which is analogous to equation \eqref{eq:Discrete_E_F(s)}.


The computations in \citep{Murphy:ADSTPF-2017,Murphy:ADSRF-2019} and those displayed in \figref{fig:BC-Flow} 
show that the spectral origin $\lambda=0$ for advection diffusion plays the role of the spectral endpoints 
$\lambda=0,1$ for two-component and polycrystalline composite materials. 
An increase in spectral mass at $\lambda=0$ increases the advection-driven 
enhancement of effective diffusivity above the 
% thermal conductivity of sea ice with no brine flux.
bare molecular diffusivity $\varepsilon$ in the advection dominated regime, 
where
$\varepsilon\ll1$ (or P\'eclet number $\xi\gg1$).
For example, the closed streamlines shown in the leftmost panel of 
\figref{fig:BC-Flow}(a) for BC-cell-flow, with fluid velocity 
field $v=(C\cos{y},B\cos{x})$ and $B=C=1$, transport tracers in a short 
range periodic motion, so long range transport is only possible due to 
molecular diffusion. Consequently, in the advection dominated regime, the effective diffusivity
scales as $\kappa^*\sim\varepsilon^{1/2}$
\citep{Fannjiang:1994:SIAM_JAM:333,Fannjiang:1997:1033,Murphy:ADSRF-2019},
vanishing as $\varepsilon\to0$. As shown in \figref{fig:BC-Flow}(b) and (c), and also in \citep{Murphy:ADSTPF-2017}, this is
reflected in the spectral measure $\nu$ by the lack of adequate mass near
$\lambda=0$ for the singular integrand $1/(\varepsilon^2+\lambda^2)$ to overcome the multiplicative factor of $\varepsilon$ for
$\kappa^*=\varepsilon(1+\langle|\nabla w_k|^2\rangle)$ in
\eqref{eq:Eff_Diffusivity_Sobolev}. 


On the other hand, when $B\ne C$ the 
streamlines elongate and connect to neighboring cells which gives rise to 
long range advection of tracers, even in the absence of molecular diffusion.
This is reflected in the spectral measure by a buildup of adequate mass near
$\lambda=0$ for the singular integrand $1/(\varepsilon^2+\lambda^2)$ to overcome the multiplicative factor of $\varepsilon$ for 
$\kappa^*=\varepsilon(1+\langle|\nabla w_k|^2\rangle)$ in
\eqref{eq:Eff_Diffusivity_Sobolev}, leading to a non-zero value of
$\kappa^*$ in the limit $\varepsilon\to0$. 
This is a key example of how the behavior of the spectral measure $\nu$ governs the behavior of the bulk transport coefficient 
$\kappa^*$.


%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering
\includegraphics[width=\textwidth] 
{RMT_melt_pond_spacings_brine_long_range.eps}
\caption{\emph{Eigenvalue spacing statistics for the sea ice melt ponds (a) and long range eigenvalue statistics for brine structures in sea ice (b).}
(a): Eigenvalue spacing distributions (ESD)  $P(z)$ for melt ponds shown in \figref{fig:melt_ponds_SM_E} corresponding to melt water area fractions 9\%, 27\%, and 57\%. 
(b): 
Spectral statistics for  brine structures  shown in \figref{fig:brine_SM_E} corresponding to area fractions of water 12\%, 51\%, and 70\%.  
We see the transition to universal Wigner-Dyson statistics as ocean phases and brine phases become connected over the scale of the sample.}
\label{fig:ESD_melt_ponds_long_range_brine}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%


\section{Random matrix theory in sea ice physics}
% \subsection{Anderson Transition for Two Phase Composites.}\label{anderson_loc}
% \section{Anderson Transition for Two Phase Composites}
% \vspace{-1ex}
%
In 
% RMT
random matrix theory (RMT)
\citep{Guhr:1998:PR.299.4.189,Bohigas:Giannoni:MR769113,Deift:2009:RMT},
long and short range 
correlations of the bulk eigenvalues away from the spectral edge 
\citep{Canali:1996:PRB.53.7.3713,Guhr:1998:PR.299.4.189}
%of matrices with random entries are
for random matrices are
measured using various eigenvalue
statistics \citep{Guhr:1998:PR.299.4.189,Bohigas:Giannoni:MR769113},
such as the 
eigenvalue spacing distribution (ESD),
spectral rigidity $\Delta_3$
and number variance $\Sigma^2$. 
To observe statistical fluctuations of these bulk eigenvalues about the
mean density, the eigenvalues must be unfolded
\citep{Bohigas:Giannoni:MR769113,Guhr:1998:PR.299.4.189,Canali:1996:PRB.53.7.3713,Plerou:2002:PRE.65.6.066126}.
The localization properties of the eigenvectors are measured in terms of quantities such as the inverse participation ratio (IPR) 
% and the eigenvector entropy 
\citep{Plerou:2002:PRE.65.6.066126,Evers:RMP:2008}. 


\begin{figure}
\centering
\includegraphics[width=\textwidth] 
{brine_anderson_localization.png}
\caption{
\emph{Anderson localization transition for electric fields in sea ice brine microstructure.}
 Examples of (a) localized and (b) extended electric fields for sea ice brine microstructure with increasing connectedness from left to right. 
The values of $s$ are taken to satisfy $\text{Im}(s) = 0.001$ with 
$0 \le \text{Re}(s) \le 1$. The electric fields are displayed for values of 
$\text{Re}(s)$ which (a) maximize and (b) minimize IPR$[E](s)$, which correspond to 
the values of such $s$ associated with (a) most localized and (b) most extended electric fields. (c) Corresponding IPR for eigenvectors $u_j$ plotted versus eigenvalue index $j$. The vertical lines define the $\delta$-components of $\mu$, where the eigenvalues satisfy $\lambda_i\lesssim10^{-14}$ and $1-\lambda_i\lesssim10^{-14}$.
The horizontal lines mark the IPR value $I_{GOE}=3/N_1$ for the Gaussian orthogonal ensemble (GOE) with matrix size $N_1\approx \phi N$, where $N=\mathrm{L}^dd$. }


%The IPR's of eigenvectors $u_j$ associated with the melt pond microstructure shown above plotted versus eigenvalue index $j$ and increasing connectedness from left to right. The vertical lines define the $\delta$-components of $\mu$
%while the horizontal lines mark the IPR value $I_{GOE}=3/N_1$ for the Gaussian orthogonal ensemble (GOE) with matrix size $N_1\approx \phi N$, where $N=\mathrm{L}^dd$.
\label{IPR_brine}
\end{figure}


In \citep{Murphy:PRL:2017}, 
we found that 
% an Anderson-like transition to universality 
% arises in the 
% %resonances of the 
% spectral characteristics of the random matrix $M=\chi_1\Gamma\chi_1$ which governs the
% transport properties of two phase media.
% %statistical properties of 
% %resonances related to the AC dielectric
% %response of percolating inhomogeneous materials. 
% %Our results establish connections between the physics of
% %Anderson localization, the mathematics of random matrix theory, 
% %and critical behavior of transport in composite materials. 
as a percolation threshold is approached and long range
order develops,
the behavior of the ESD transitions from weakly-correlated Poissonian statistics 
toward obeying universal Wigner-Dyson (WD) statistics of the Gaussian Orthogonal
Ensemble (GOE).
% , as illustrated in \figref{fig:ESD} %\figref{spectral_measure}.
The eigenvectors de-localize, and
mobility edges appear \citep{Murphy:PRL:2017}, similar to the metal/insulator transition in 
solid state physics.
% , as shown in \figref{Eigenvector_Localization_RRN}.
%statistical properties of 
%resonances related to the AC dielectric
%response of percolating inhomogeneous materials. 
%Our results establish connections between the physics of
%Anderson localization, the mathematics of random matrix theory, 
%and critical behavior of transport in composite materials. 
We explored the transition 
in the 2D and 3D 
RRN,
% random resistor network (RRN), 
as well as in
sea ice microstructures such as
%that associated with the effective complex permittivity $\epsilon^*$
%for  
in 2D discretizations of the brine microstructure of sea ice 
\citep{Golden:S-2238,Golden:GRL:L16501,Golden:NAMS:2009}, 
melt ponds on Arctic sea ice
\citep{Hohenegger:TC-2012}, 
% as shown in \figref{spectral_measure} (d)--(f), 
the sea ice pack itself, and 
porous 
human
bone \citep{Golden:JBM:337,Kabel:B:1999,Bonifasi-Lista:PMB-3063, Cherkaev:JBiomech-345}. 
We extended these results to two-component composite media with quasiperiodic microgeometry in \citep{Morison:2022:quasiperiodic}.











%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Description of how we created the matrix $G=\chi_1\Gamma\chi_1$ %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To calculate the spectral measures for 2D images of sea ice and bone microstructure, we converted photographs to binary images corresponding to the two phases. For example, \figref{spectral_measure} (d)--(f) displays binary images of Arctic melt ponds with varying degrees of connectivity. 
% The 
% % binary 
% images were then mapped to high density resistor networks, 
% from which the matrix 
% $M$
% % $M=\chi_1\Gamma\chi_1$ 
% was then obtained.
% %for these discretizations.






% MUCH MORE DETAILS IN THIS PARAGRAPH THAN THE ABOVE PARAGRAPH !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% Computing the spectral measure $\mu$ for a given composite microstructure first involves discretizing a two phase image of the composite into a square lattice filled with 1's and 0's corresponding to the two phases. Then, we can compute the spectral measure of the matrix $M=\chi_1\Gamma\chi_1$ -- depending on the geometry via $\chi_1$ -- directly from the eigenvalues $\lambda_i$ and eigenvectors $\vecw_i$ of this matrix. 
% X-ray CT renderings of brine microstructure and aerial photographs of Arctic melt ponds and sea ice floes, as well as human bone microstructure were converted to 2D binary fluid/ice representations which describe the geometry and the component connectivity of the system. Statistical realizations of these binary composite structures determine the statistical realizations of the associated random matrix. The statistical properties of the computed eigenvalues and eigenvectors of the random matrix provide valuable information regarding the behavior of the associated spectral measure. 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Detailed description of Anderson transition for RRN %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The measure $\mu$ exhibits fascinating transitional behavior as a
% function of system connectivity.
% % , which controls critical behavior of $\epsilon^*$ near connectedness thresholds.
% For example, in the case
% of a 
% RRN
% %random resistor network (RRN) 
% with a low 
% volume fraction $p$ of open bonds, as shown in
% \figref{spectral_measure}a,
% %Fig.\,1A,
% there are spectrum-free
% regions at the spectral endpoints
% $\lambda=0,1$~\citep{Murphy:JMP:063506}. However, as  
% $p$ approaches the 
% percolation threshold $p_c$~\citep{Stauffer-92,Torquato:RHM-02}
% and the system becomes increasingly connected, these spectral gaps shrink and
% then vanish~\citep{Murphy:JMP:063506,Jonckheere_Luck_JPA_1998}, 
% as shown in \figref{spectral_measure}a--c, 
% leading
% to the formation of $\delta$-components of $\mu$ at the spectral endpoints,
% \emph{precisely}~\citep{Murphy:JMP:063506} when $p=p_c$ (and
% $p=1-p_c$ in $d=3$). This leads to critical behavior of $\epsilon^*$ for
% insulating/conducting and conducting/superconducting
% systems \citep{Murphy:JMP:063506}. 
% This gap behavior of $\mu$ has 
% led \citep{Golden:PRL-3935,Murphy:JMP:063506}
% %~\citep{Golden:PRL-3935,Murphy:JMP:063506}
% to a detailed description
% of these critical transitions in $\epsilon^*,$ which is directly analogous to
% the Lee--Yang--Ruelle--Baker
% description \citep{Baker-1990,Golden:PRL-3935} of the 
% Ising model phase transition in the magnetization.
% % $M$. 
% Moreover, using this gap behavior, all of the classical critical
% exponent scaling relations
% %~\citep{Bergman:SSP-147}
% were recovered~\citep{Murphy:JMP:063506,Golden:PRL-3935} without heuristic
% scaling forms but instead by using 
% the \emph{rigorous} integral representation for
% $\epsilon^*$ in \eqref{eq:F(s)} involving $\mu$. 
% %This spectral behavior emerges in all the systems mentioned above.
% % variety of composite systems, such as the brine microstructure of sea ice~\citep{Golden:S-2238,Golden:GRL:L16501,Golden:NAMS:2009}, melt
% % ponds on the surface of Arctic sea ice~\citep{Hohenegger:TC-2012},
% % the sea ice pack itself
% % (as shown in \figref{spectral_measure}a--c),
% % and porous human bone~\citep{Kabel:B:1999}.



% We show that this spectral behavior emerges in a variety of composite
% systems, such as the brine microstructure of sea
% ice~\citep{Golden:S-2238,Golden:GRL:L16501,Golden:NAMS:2009}, melt
% ponds on the surface of Arctic sea ice~\citep{Hohenegger:TC-2012}, the
% sea ice pack itself, and porous human
% bone~\citep{Kabel:B:1999}. Our results indicate that it is pervasive in
% such macroscopic systems 
% %They span a broad range of application areas as well as length scales,
% %and
% and arises simply from connectedness -- at
% the most basic level of characterizing any physical system with
% inhomogeneities. 



% \begin{figure}%[h!] 
% \centering
% \includegraphics[width=\textwidth]{Eigenvector_Localization_RRN_Electric_Fields_RRN_Brine_Bone.jpg}
% \caption{\textbf{Spectral behavior of composite structures.}
% %   (\textbf{a}) 
% %   %Displayed from left to right are the full electric field $\chi_1 \vecE$ (in log scale) and 
% %   Examples of localized and extended eigen-modes of
% %   $\chi_1\vecE$ 
% %   %(in linear scale) 
% %   for a realization of the 2D
% %   random resistor network (RRN), with a system size $\mathrm{L}=50$,
% %   volume fraction $p=p_c=1/2$, 
% %   and $\epsilon_1/\epsilon_2=4.0\times10^{10}$.    
% %   , and the DC values of the component conductivities are those for 
% %   silver and silicon at 20 $^\circ$C, $\epsilon_1=6.3\times10^7$ S/m and 
% %   $\epsilon_2=1.56\times10^{-3}$ S/m. 
% %   The corresponding localized and extended eigenvectors have inverse 
% %   participation ratios (IPR's) $I_n\approx0.081$ and $I_n\approx0.002$
%   (a) The ESDs for 3D RRN with $p_c\le p\le 1-p_c$.
%   (b) The IPR's of eigenvectors $u_j$ associated with realizations of the
%   3D RRN plotted versus eigenvalue index $j$, with $\mathrm{L}=12$
%   and increasing values of $p$ from left to right, with
%   $1-p_c\approx0.7512$. The vertical lines define the $\delta$-components of $\mu$
% %   spectral measure $\mu$ at the left and right spectral endpoints, where
% %   the associated eigenvalues $\lambda_n$ satisfy $\lambda_n\lesssim10^{-14}$ and
% %   $\lambda_n\gtrsim1-10^{-14}$, respectively, 
%   while the horizontal lines mark the 
%   IPR value 
%   $I_{GOE}=3/N_1$ 
%   for the Gaussian orthogonal ensemble (GOE)
%   with matrix size $N_1\approx pN$, where $N=\mathrm{L}^dd$.
%   (c) The
%   ensemble averaged IPR $\langle I\rangle$ as a function of $p$, displaying
%   transitional behavior at the percolation thresholds, $p_c=1/2$ for 2D
%   and $p_c\approx0.2488$ for 3D.
% %   (\textbf{d}) The normalized, ensemble averaged
% %   entropy $\langle S\rangle/S_{GOE}$ and associated localization length
% %   $\ell/\ell_{GOE}$ as a function of $p$, where $S_{GOE}\approx\ln(N_1/2.07)$ and
% %   $\ell_{GOE}=N_1$.
% (d)--(f) Electric fields (with logarithmic scale) in composite structures calculated via \eqref{eq:Electric_Field_two_component}.
% %   Panels (b)--(d) demonstrate that, as $p$ increases
% %   and the systems become increasingly connected, the eigenvectors
% %   become progressively extended, on average, with decreasing 
% %   $\langle  I\rangle$.
%   %and increasing $\langle S\rangle/S_{GOE}$ and $\ell/\ell_{GOE}$.       
% }
% \label{Eigenvector_Localization_RRN}      
% \end{figure}
%




% Our computations of the ESD in \figref{spectral_measure} reveal a mechanism for the collapse in the spectral gaps of $\mu$. 
For highly correlated WD spectra 
exhibited by, for example, real-symmetric matrices of the 
%Gaussian orthogonal ensemble (GOE),
GOE,
the nearest neighbor 
%eigenvalue spacing distribution (ESD) $P(z)$
ESD $P(z)$
is accurately
approximated by $P(z)\approx(\pi z/2)\exp(-\pi z^2/2)$ 
%(Wigner's surmise),
which illustrates \emph{eigenvalue repulsion} and
vanishes linearly as spacings $z \rightarrow 0$ 
\citep{Guhr:1998:PR.299.4.189,Stone:1991:chapter9:369,Canali:1996:PRB.53.7.3713}. 
In contrast, the ESD for uncorrelated Poisson spectra, 
$P(z)=\exp(-z)$,
allows for significant level degeneracy \citep{Guhr:1998:PR.299.4.189}. 
In \figref{fig:ESD_melt_ponds_long_range_brine}(a) we display the ESDs 
for Poisson (blue) and WD (green) spectra, along with the behavior of 
the ESDs for the matrix $M=\chi_1\Gamma\chi_1$, 
corresponding to the Arctic melt ponds in
\figref{fig:melt_ponds_SM_E} with fluid area fraction $\phi$. 
It shows
that for sparsely connected systems, the behavior of the ESDs is
well described by weakly correlated
% Poisson-like
Poisson-like
statistics \citep{Canali:1996:PRB.53.7.3713}. 
% They increase linearly from zero but the initial 
% slope of the curve is steeper than in WD case, implying less repulsion, and the tails decay exponentially. 
With increasing connectedness, 
%and long range order is established, 
the ESDs transition toward highly correlated WD statistics with strong
level repulsion. 
%and Gaussian tails.
% The blue dash-dot curve displayed in \figref{fig:Spacings}
% is the ESD for Poisson spectra, while the green
% dashed curve is the ESD for the GOE.
This behavior of the ESD reveals a mechanism for the collapse in the
spectral gaps of $\mu$. For sparsely connected systems, the weak level
repulsion allows for significant level degeneracy and resonances in $\mu$
as shown in \citep{Murphy:2015:CMS:13:4:825} for the 2D percolation model
and in \figref{fig:melt_ponds_SM_E} (a) for Arctic melt pond microstructure. As
the system becomes increasingly connected, the level repulsion
increases, causing the eigenvalues to spread out, which, in turn, causes the
gaps in the measure near the spectral edges to collapse and subsequently 
form $\delta$-components of the measure at the spectral endpoints 
$\lambda=0,1$ (as shown in the figures in \secref{sec:two_component_discrete}).
Our computations
of $\Delta_3$ and $\Sigma^2$ are are shown in
\figref{fig:ESD_melt_ponds_long_range_brine}(b) for the brine
microstructure in \figref{fig:brine_SM_E}. Notice that as the system becomes increasingly connected, it transitions towards that of the GOE, indicating an 
increase in the long range correlations of the eigenvalues.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Description of the symmetry in the eigenvalue density and ESD %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For the 2D and 3D RRN, the eigenvalue
% density $\rho(\lambda,p)$ obeys 
% $\rho(\lambda,p)=\rho(1-\lambda,1-p)$ in the
% \emph{bulk} of the spectrum. This is reflected in the ESDs by the
% symmetry $P(z,p)=P(z,1-p)$, as shown for the 2D RRN in
% \figref{spectral_measure}a-c. 
% Our computations of 
% % the spectral rigidity
% $\Delta_3$ 
% % and
% % and number variance
% % $\Sigma^2$
% %which probe long range eigenvalue correltions, 
% exhibit
% % exhibit
% a similar transitional behavior \citep{Murphy:PRL:2017}.
% The ESD's displayed in \figref{fig:Spacings}d for the 3-dimensional
% percolation model suggest that the GOE limit is attained for all $p_c\leq
% p\leq1-p_c$. 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Relation to Anderson MIT in 3D %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Our computations of the ESD, 
% $\Delta_3$ and $\Sigma^2$
% % spectral rigidity, and number variance 
% for 3D RRN also suggest that the GOE limit is attained 
% by these eigenvalue statistics for
% %by the long  range statistics for the 3D RRN for 
% \emph{all} $p_c\leq p\leq1-p_c$, appearing to overlie the GOE limit almost exactly for all  $p$ values tested in $p_c \leq p \leq 1-p_c$ \citep{Murphy:PRL:2017}.
% % as shown for ESDs in \figref{Eigenvector_Localization_RRN}a.
% % The computed ESDs also 
% % appear to overlie the GOE limit almost exactly for all
% % $p$ values tested in $p_c \leq p \leq 1-p_c$.
% %The 2D RRN displays a similar behavior for $p=p_c=1/2$.  
% With this in mind, we recall the  
% Anderson transition, where low disorder corresponds to 
% extended states and WD statistics. When disorder exceeds a critical level,
% the states localize and the eigenvalues become de-correlated.
% We view the 3D RRN with $p_c \leq p \leq 1-p_c$ to be ``ordered''
% with extended states and WD statistics. As $p$ decreases, the
% disorder $-$ or blockages to the flow $-$ increases, and the eigenstates 
% localize. % as we now discuss.
% % and the eigenvalue repulsion diminishes.
% %Hence for transport in two-phase composites, $p$ plays the role of order %parameter.
% %This fascinating spectral behavior did not emerge discovered in related numerical investigations \citep{Jonckheere_Luck_JPA_1998,
	% %Gu:2002:PRE.65.046129,Gu:2002:PRB.66.012202}


The eigenvectors $u_j$ of $M=\chi_1\Gamma\chi_1$,
% $j=1,\ldots,N_1$, 
associated with the %random matrix 
$N_1\times N_1$ sub-matrices of $\Gamma$, 
%discussed above, 
also exhibit a connectedness driven transition in their localization properties. 
The IPR is %a commonly used quantity to describe the 
%localization of eigenvectors \citep{Plerou:2002:PRE.65.6.066126}.
% The IPR
% $I_j$ \citep{Plerou:2002:PRE.65.6.066126} is 
% %a commonly used quantity
% %to measure the localization of
% % the eigenvector
% %$\vecv_j$. It is
defined as $I_j=\sum_i|u_j^i|^4$, $i,j=1,\ldots,N_1$, where $u_j^i$ is the $i$th component of $u_j$. 
%Two limiting cases 
% illustrate the meaning of $I_j$:
%  $(i)$ a normalized vector with only one component
% $v_j^i=1$ has $I_j=1$; $(ii)$ a vector with identical
% components $v_j^i=1/\sqrt{N_1}$ has $I_j=1/N_1$. 
Eigenvectors of
matrices in the GOE are known to be highly extended 
% and 
%uniformly distributed on the unit sphere in
%$\mathbb{R}^N$~\citep{Guhr:1998:PR.299.4.189} with statistics  
%characterized by the squared Gaussian $\chi^2$ distribution
%law~\citep{Porter:1965:Book:Stat_Fluct}, 
% independent of the eigenvalue distribution 
%of the eigenvalues
\citep{Deift:2009:RMT}, 
%and having Gaussian distributed
%amplitudes~\citep{Bermudez:2015:PRE.91.032122,Plerou:2002:PRE.65.6.066126}.
%In this case,
with asymptotic value of  
the IPR given by 
$I_{GOE}=3/N_1$ \citep{Plerou:2002:PRE.65.6.066126}.
% 
% with large magnitudes in only a few
% %type-one
% resistors, while extended eigenvectors correspond
% to fields extending throughout the
% network.
% \figref{Eigenvector_Localization_RRN}d-f shows our computations of the electric field $\chi_1 \vecE$ for various composite structures.
%he 2D RRN with $p=p_c=1/2$. 
%
%
%  \figref{Eigenvector_Localization_RRN}b 
% displays $I_j$
% %$j=1,\ldots,N_1$, 
% for the 3D RRN 
% %$p=1-p_c\approx0.7512$,
% %$p=1-p_c$,
% plotted \textit{vs.} index $j$, for various $p$ values.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Details about the delta functions forming at p_c %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% In numerical simulations of finite RRN the $\delta$-component at  
% %the left spectral endpoint
% $\lambda=0$ is manifested 
% by a large number of $\lambda_j$ with magnitude
% $\lesssim10^{-14}$, followed by an abrupt 
% change of magnitude $\gtrsim10^{-4}$,
% with no eigenvalues in the interval $(10^{-14},10^{-4})$, 
% and similarly for
% %the right spectral endpoint
% $\lambda=1$.
% % The locations of abrupt changes in
% % eigenvalue magnitudes are identified by red \emph{vertical} lines.
% % The GOE IPR value is identified by the red
% % \emph{horizontal} line. This %(and others like it) 
% % demonstrates
% % that eigenvectors with the $\delta$-components
% % at $\lambda=0,1$
% % are %typically 
% % more extended than others, 
% % with $I_j$ values  
% % closer to the GOE limit. 
% %
In \citep{Murphy:PRL:2017}, 
we found for the 2D and 3D percolation models that as $p$ surpasses $p_c$
and long range order is established in a RRN, \textbf{``mobility edges''}
form in the eigenvector IPR with a sudden increase in the number of
extended eigenvectors. 
%This behavior in the random matrix
%$G=\chi_1\Gamma\chi_1$ 
This is analogous to Anderson localization, 
where mobility edges mark the characteristic energies of the metal-insulator transition 
(MIT) \citep{Guhr:1998:PR.299.4.189}.
Remarkably, the mobility edges for RRN
% in \figref{Eigenvector_Localization_RRN}b
are due to very extended eigenstates associated with $\delta$-components
%at $\lambda=0,1$,
% (red \emph{vertical} lines) 
that form at the spectral endpoints \emph{precisely} at the percolation threshold $p_c$ (and $1-p_c$ for 3D) \citep{Murphy:JMP:063506},
which control critical behavior in 
insulator/conductor and conductor/superconductor
systems \citep{Murphy:JMP:063506,Clerc:AP-191,Bergman:SSP-147}.  This and other eigenvector phenomena was observed for two-component composite media with quasiperiodic microgeometry in \citep{Morison:2022:quasiperiodic}.



The IPR phenomena for the eigenvectors of the matrix $G=\chi_1\Gamma\chi_1$ is shown for sea ice brine microstructure in \figref{IPR_brine} (c). The electric field $E$ within sea ice brine microstructure exhibits a frequency dependent Anderson localization transition, as shown in \figref{IPR_brine} (a) and (b). To generate these figures, the IPR of the electric field, $IPR[E](s)$ in the brine phase $\chi_1 E$ was calculated as a function of $s$ via equation \eqref{eq:Discrete_E_F(s)} for $\text{Im}\,s=0.001$ and $0\le\text{Re}\,s\le1$, then normalized to have unit length. The values of $s$ were then selected where $IPR[E](s)$ attains its maximum and minimum, corresponding to the most localized and most extended $E$, respectively. For those values of $s$, the electric fields for brine microstructure are shown in \figref{IPR_brine} (a) and (b). The localized electric fields in (a) are characterized by high intensity "hot spots" at the brine-ice interfaces, while extended electric fields have high intensities spread out more evenly across the connected and near-connected brine components. 




% \begin{figure}[t] \centering
	% \includegraphics[width=6.3in]{anderson_localization_pack_ice_connected_L_200_log.jpg}
	% \caption{\textbf{Anderson localization of horizontal electromagnetic waves through a composite geometry from pack ice.} (a) Spectral function $\mu(\lambda)$ (histogram representation of the spectral measure $\mu$) plotted \textit{vs.} the spectral variable $\lambda$. The area under the curve is the area fraction of the ocean phase. The strong resonance at $\lambda=0$ indicates the ocean phase is connected on large scales. (b) The inverse participation ratio (IPR) for the eigenvectors $\vecw_i$ of the matrix $G$ plotted \textit{vs.} index $i$. The $1/N$ and $3/N$ IPR values for completely extended and GOE vectors are shown as horizontal black and green lines.
		% %, respectively. 
		% The red vertical lines indicate Dirac-$\delta$ components of $\mu$ at $\lambda=0,1$, numerically manifested as a large number of isolated eigenvalues satisfying $\lambda_i\lesssim10^{-13}$ and $1-\lambda_i\lesssim10^{-13}$, %respectively, 
		% and mark eigenvector ``mobility edges.'' (c) The IPR for the point-wise modulus of the electric field $|\vecE|$ (equivalently current $|\vecD|$) normalized to have unit length and plotted \textit{vs.} $Re\,s$ (equivalently frequency $\omega$) for $Im\,s=-0.001$. Larger and smaller values of $IPR(|\vecE|)$ indicate frequencies where $|\vecE|$ is more localized and extended, respectively, being especially extended for $\lambda=0,1$. (d)--(e) For the same floe geometry, some frequencies of the applied electromagnetic wave lead to primarily localized fields within the ice pack (d) or fields extended throughout the pack (e), shown in log scale to reveal detail.
		% }  
	% \label{fig:anderson_localization}
	% \end{figure}
% %




\conclusions  %% \conclusions[modified heading if necessary]
We have given a tour through various problems of sea ice 
physics concerned with homogenization and 
how they can be rigorously addressed with 
the powerful analytic continuation method and its extensions. 
The effective complex
permittivity of sea ice treated as a two-phase composite of
pure ice with brine inclusions, or treated as a polycrystalline
material, and the effective diffusivity for advection diffusion
problems are all Stieltjes functions of their variables.
We showed how these functions have integral representations 
involving spectral measures which distill the mixture or
velocity field geometries into the spectral properties of
a self adjoint operator like the Hamiltonian in quantum physics.
These spectral representations have been used to obtain 
rigorous forward and inverse bounds
on effective transport coefficients for sea ice. Finally, the combination of random matrix theory with these spectral representations has revealed fascinating parallels between quantum transport in 
disordered media, Anderson localization, and classical transport phenomena.

%% The following commands are for the statements about the availability of data sets and/or software code corresponding to the manuscript.
%% It is strongly recommended to make use of these sections in case data sets and/or software code have been part of your research the article is based on.

%\codeavailability{TEXT} %% use this section when having only software code available


%\dataavailability{TEXT} %% use this section when having only data sets available


%\codedataavailability{TEXT} %% use this section when having data sets and software code available


%\sampleavailability{TEXT} %% use this section when having geoscientific samples available


%\videosupplement{TEXT} %% use this section when having video supplements available


\appendix
% \section{}    %% Appendix A

% \subsection{}     %% Appendix A1, A2, etc.


% \section{Existence of field decompositions}\label{sec:field_decomp}
% %
% In this section, for the discrete setting in Sections \ref{sec:two_component_discrete} 
% and \ref{polycrystal_discrete}, we prove that there exists an electric field 
% $\vecE$ satisfying discrete versions of equations 
% \eqref{eq:Maxwells_Equations} and \eqref{eq:E_variational}. Towards this goal,
% we follow \citep{Huang:2019:AMSA.4.1} and consider the finite difference 
% representations of the partial differential operators $\partial_i\to C_i$,
% $i=1\ldots,d$, where $d$ denotes dimension. The matrices $C_i$ depend on boundary
% conditions which, without loss of generality, we take to be periodic boundary 
% conditions. Denote the matrix representation 
% of the gradient operator (using Matlab vertical block notation) 
% by $\nabla=[C_1;\ldots;C_d]$. The discretization of the divergence operator is 
% given by $-\nabla^T$ and the discrete curl operator is given by
% \citep{Huang:2019:AMSA.4.1}
% %
% \begin{align}\label{eq:block_Sigma}
%   C&=\left[
%   \begin{array}{rrr}
%     O&-C_3&C_2\\
%     C_3&O&-C_1\\  
%     -C_2&C_1&O
%     \end{array}
%     \right]
% \;\text{ for 3D,} 
% \\
% C&=[-C_2,C_1]
% \;\text{ for 2D.} 
% \notag
% %
% \end{align}
% %
% The operators $C_i$, $i=1,2,3$, in \eqref{eq:block_Sigma} are normal 
% and commute with each other \citep{Huang:2019:AMSA.4.1},
% %
% \begin{align}
%     C_i^T C_j=C_jC_i^T \text{ and } C_iC_j=C_jC_i, \text{ for } i,j=1,2,3.
% \end{align}
% %
% Consequently, the discrete form of equations \eqref{eq:Maxwells_Equations} 
% and \eqref{eq:E_variational}, which we'll establish in this section, is
% %
% \begin{align}\label{eq:Maxwells_Equations_discrete_decomp}
% C\vecE=0, 
% \quad -\nabla^TJ=0, 
% \quad J=\beps E,
% \quad
% \vecE=\vecE_f+\vecE_0\,,
% \quad
% C\vecE_f=0\,,    
% \quad
% \langle\vecJ\cdot\vecE_f\rangle=0\,,
% \quad
% \langle\vecE_f\rangle=0\,,
% \end{align}
% %
% where in this finite size discrete setting, 
% $\langle\cdot\rangle$ denotes volume average followed by ensemble average
% \citep{Murphy:2015:CMS:13:4:825,Murphy:JMP:2020}.
% To set notation, denote by $\text{Ran}(B)$ and $\text{Ker}(B)$ 
% the range and kernel (null space) of the matrix $B$, respectively 
% \citep{Horn_Johnson-1990}.  
% Therefore, we seek to prove there exists a vector $E$ 
% satisfying $E\in\text{Ker}(C)$ such that $E=E_f+E_0$, where 
% $E_f\in\text{Ker}(C)$ and $\langle E_f\rangle=0$ so 
% $\langle\vecE\rangle=\vecE_0$. 
% Moreover, we seek to find a vector 
% $\vecJ\in\text{Ker}(\nabla^T)$ satisfying $\vecJ=\beps\vecE$ and 
% $\langle\vecJ\cdot\vecE_f\rangle=0$.


% We now summarize some useful identities relating the discrete representations
% of the gradient, divergence, and curl operators which follow from these properties 
% of the matrices $C_i$ \citep{Huang:2019:AMSA.4.1},
% %
% \begin{align}\label{identities}
% &\Delta=\nabla\cdot\nabla\to-\nabla^T\nabla\,,
% \\
% &\bDelta = \text{diag}(\Delta,\ldots,\Delta)
%        \to I_d\otimes(\nabla^T\nabla)\,,
% \notag\\
% &\nabla\times\nabla\times\to C^TC\,,
% \notag\\
% &\nabla\times\nabla\times=\nabla(\nabla\cdot)-\bDelta
% \to -\nabla\nabla^T + I_d\otimes(\nabla^T\nabla)\,,
% \notag\\
% &\nabla\cdot(\nabla\times)\to -\nabla^TC^T=0\,,
% \notag\\
% &\nabla\times\nabla\to C\nabla=0\,,
% \notag
% \end{align}
% %
% where $\otimes$ denotes the Kronecker product.
% The last two identities $\nabla^TC^T=0$ 
% and $C\nabla=0$ in equation \eqref{identities} indicate that 
% % $\text{Ran}(C^T)\subseteq\text{Ker}(\nabla^T)$ and 
% % $\text{Ran}(\nabla)\subseteq\text{Ker}(C)$. 
% %
% \begin{align}\label{subspaces}
% \text{Ran}(C^T)\subseteq\text{Ker}(\nabla^T)\,,
% \quad
% \text{Ran}(\nabla)\subseteq\text{Ker}(C)
% \end{align}
% %
% We can now restate our goal in \eqref{eq:Maxwells_Equations_discrete_decomp} as, 
% find "potentials" $\varphi$ and $\psi$ such that the vectors 
% $\vecE$ and $\vecJ$ in equation 
% \eqref{eq:Maxwells_Equations_discrete_decomp} satisfy $E=E_f+E_0$ with 
% $E_f=\nabla\varphi$ and $E_0\in\text{Ker}(\nabla)$, and $J=C^T\psi+J_0$, 
% where $J_0\in\text{Ker}(C^T)$.
% The last two identities \eqref{identities} provide a relationship between 
% rank and kernel of the operators $C$, $\nabla$, and their transposes. 
% The fundamental theorem of linear algebra provides a relationship 
% between the range of a matrix $A$ and the kernel of it's transpose $A^T$, 
% which will be useful later in this section.
% %
% \begin{theorem}[Fundamental theorem of linear algebra]\label{FTLA}
% %
% Let $A$ be a real valued matrix of size $m\times n$ then
% %
% \begin{align} \label{eq:FTLA}
% \mathbb{R}^m=\text{\emph{Ran}}(A)\oplus\text{\emph{Ker}}{(A^T)}\,,
% \qquad
% \mathbb{R}^n=\text{\emph{Ran}}(A^T)\oplus\text{\emph{Ker}}{(A)}\,,
% \end{align}
% %
% where $\oplus$ indicates $\text{Ran}(A)$ is orthogonal 
% to $\text{Ker}{(A^T)}$, i.e., $\text{Ran}(A)\perp\text{Ker}{(A^T)}$, 
% for example.
% \end{theorem}
% %




% Applying Theorem \ref{FTLA} to the matrices $\nabla$ 
% and $C^T$ indicates that
% $\mathbb{R}^m=\text{Ran}(\nabla)\oplus\text{Ker}{(\nabla^T)}$ and 
% $\mathbb{R}^m=\text{Ran}(C^T)\oplus\text{Ker}{(C)}$. Therefore,
% from equation \eqref{subspaces} we have that
% divergence-free fields are orthogonal to gradients (curl-free fields)
% and curl-free fields are orthogonal to $\text{Ran}(C^T)$ (divergence-free fields). This is a discrete version of the Helmholz Theorem, which states that curl-free and divergence-free fields 
% (or, in other words, the gradient and cycle spaces) are mutually 
% orthogonal. This also establishes the important relationship 
% % $\text{Ran}(C^T)\perp\text{Ran}(\nabla)$. 
% %
% \begin{align}\label{eq:perp_C_grad}
% \text{Ran}(C^T)\perp\text{Ran}(\nabla)\,.    
% \end{align}
% %





% Orthogonal bases can be given for each of 
% the mutually orthogonal spaces in equation \eqref{eq:FTLA} through the singular
% value decomposition (SVD) \citep{Horn_Johnson-1990} of the matrix $A=U\Sigma V^T$,
% which also provides important information relating the matrices $C$, $\nabla$, etc. 
% Here $U$ and $V$ are 
% orthogonal matrices of size $m\times m$ and $n\times n$ satisfying $U^T U=UU^T=I_m$ 
% and $V^T V=VV^T=I_n$, where $I_m$ is the identity matrix of size $m$. Moreover, $\Sigma$
% is a diagonal matrix of size $m\times n$ with diagonal components consisting of the 
% positive \emph{singular values} $\nu_i$, $i=1,\ldots,n$, of the matrix $A$ satisfying 
% $\nu_1\ge\nu_2\ge\cdots\ge\nu_\rho>0$ and $\nu_{\rho+1}=\nu_{\rho+2}=\cdots=\nu_n=0$,
% where $\rho$ is the rank of $A$. Writing the matrix $\Sigma$ in block form we have
% %
% \begin{align}\label{eq:block_Sigma}
%   \Sigma=\left[
%   \begin{array}{cc}
%     \Sigma_1&O_1\\
%     O_1^T&O_2\\  
%     O_3&O_4
%     \end{array}
% \right]\,,
% \end{align}
% %
% where $\Sigma_1$ is a diagonal matrix of size $\rho\times\rho$ with diagonal
% consisting of the strictly positive singular values, $O_1$ and $O_2$, are matrices 
% of zeros of size $\rho\times(n-\rho)$ and $(n-\rho)\times(n-\rho)$, and the bottom block 
% of zeros $[O_3,O_4]$ is of size $(m-n)\times n$. Write the matrices $U$ and $V$ in block form as
% $U=[U_1,U_0,U_2]$ and $V=[V_1,V_0]$, where $U_1$ and $V_1$ are the 
% columns of $U$ and $V$ corresponding to the strictly positive singular values in $\Sigma_1$,
% $U_0$ and $V_0$ correspond to the the singular values with value zero, $\nu_i=0$, and $U_2$
% corresponds to the bottom block of zeros $[O_3,O_4]$ in $\Sigma$. Taking in account all the
% blocks of zeros in $\Sigma$, we can write $A=U_1\Sigma_1V_1^T$. 
% We now state a well known fact about the SVD of the matrix $A$ \citep{Horn_Johnson-1990}.
% %
% \begin{align}\label{svd_col_space}
% \text{Range}(A)=\text{Col}(U_1),
% \quad
% \text{Ker}(A)=\text{Col}(V_0),
% \quad
% \text{Range}(A^T)=\text{Col}(V_1),
% \quad
% \text{Ker}(A^T)=\text{Col}([U_0,U_2]),
% \end{align}
% %
% where $\text{Col}(B)$ denotes the column space of the matrix $B$, i.e., 
% the space spanned by the columns of $B$. 




% Applying the SVD
% to the matrices $\nabla=U^\times\Sigma^\times[V^\times]^T$ 
% and $C^T=U^\bullet\Sigma^\bullet[V^\bullet]^T$ and using the orthogonality 
% of the columns of the matrices $U_1$ and $V_1$ and 
% the invertability of $\Sigma_1$, from $C\nabla=0$ in 
% \eqref{identities} we have $[U^\bullet]^T U_\times=0$, and similarly 
% $\nabla^TC^T=0$ implies $[U^\times]^T U_\bullet=0$. The formula
% $[U^\bullet]^T U_\times=0$ is a restatement of equation 
% \eqref{eq:perp_C_grad}. Writing 
% $U^\times=[U_1^\times,U_0^\times,U_2^\times]$ and 
% $U^\bullet=[U_1^\bullet,U_0^\bullet,U_2^\bullet]$ we have 
% established the $\text{Ran}(\nabla)=\text{Col}(U_1^\times)$,
% $\text{Ran}(C^T)\subseteq\text{Ker}(\nabla^T)
% =\text{Col}([U_0^\times,U_2^\times])$. Also, since
% $\text{Ran}(C^T)=\text{Col}(U_1^\bullet)$ and 
% $\text{Ran}(C^T)\perp\text{Ran}(\nabla)$, we can 
% write
% %
% \begin{align}
% U^\times=[U_1^\times,U_0^{\times\bullet},U_1^\bullet]\,,
% \quad
% U^\bullet=[U_1^\bullet,U_0^{\times\bullet},U_1^\times]\,,
% \end{align}
% %
% where the columns of $U_0^{\times\bullet}$ are orthogonal to both
% the $\text{Ran}(C^T)$ and the $\text{Ran}(\nabla)$. Since
% $U^\times [U^\times]^T=I_m$ we have the \emph{resolution of 
% the identity}
% %
% \begin{align}\label{eq:res_of_identity}
% % I_m
% % =U_1^\times[U_1^\times]^T + 
% % U_0^{\times\bullet}[U_0^{\times\bullet}]^T + 
% % U_1^\bullet[U_1^\bullet]^T
% \Gamma_\times+\Gamma_0+\Gamma_\bullet
% =I_m,
% \qquad
% \Gamma_\times=U_1^\times[U_1^\times]^T\,, 
% \quad
% \Gamma_\bullet=U_1^\bullet[U_1^\bullet]^T\,,
% \quad
% \Gamma_0=U_0^{\times\bullet}[U_0^{\times\bullet}]^T\,,
% \end{align}
% %
% where $\Gamma_\times$, $\Gamma_\bullet$, and $\Gamma_0$, 
% are mutually orthogonal projection matrices onto $\text{Ran}(\nabla)$, 
% $\text{Ran}(C^T)$ and the orthogonal complement of 
% $\text{Ran}(\nabla)\cup\text{Ran}(C^T)$. We are now ready to 
% state the main result of this section as the following theorem.
% %
% \begin{theorem}\label{thm:E_J_decomp}
% %
% Let the electric and current fields $\vecE$ and $\vecJ$ satisfy
% %
% \begin{align}\label{eq:Maxwells_Equations_discrete}
% C\vecE=0, 
% \quad -\nabla^T J=0, 
% \quad J=\beps E.
% \end{align}
% %
% Then, there exists a "potential" $\varphi$ and vector $E_0$ such that
% $E=E_f+E_0$, where $E_f=\nabla\varphi$,
% $I_d\otimes(\nabla^T\nabla))E_0=0$, $CE_f=0$, 
% $\langle\vecJ\cdot\vecE_f\rangle=0$, and $\langle\vecE_f\rangle=0$.
% \end{theorem}
% %
% \begin{proof}
% From the resolution of the identity in equation \eqref{eq:res_of_identity}
% we have $E=(\Gamma_\times+\Gamma_0+\Gamma_\bullet)E$.
% Since $\Gamma_\bullet$ projects onto the $\text{Ran}(C^T)$,
% $\mathbb{R}^m=\text{Ran}(C^T)\oplus\text{Ker}{(C)}$, and 
% $\vecE\in\text{Ker}{(C)}$ we have $\Gamma_\bullet E=0$. Denoting
% by $E_f=\Gamma_\times E$, since 
% $U_1^\times=\nabla V_1^\times[\Sigma_1^\times]^{-1}$, we can write 
% $E_f=\nabla\varphi$, where
% $\varphi=V_1^\times[\Sigma_1^\times]^{-1}[U_1^\times]^TE$. Denote 
% $E_0=\Gamma_0E$. Since $\Gamma_\times E_0=0$, $\Gamma_\times$ is a 
% projection onto $\text{Ran}(\nabla)$, and 
% $\mathbb{R}^m=\text{Ran}(\nabla)\oplus\text{Ker}{(\nabla^T)}$, we have 
% $E_0\in\text{Ker}(\nabla^T)$. Similarly, since $\Gamma_\bullet E_0=0$,
% $\Gamma_\bullet$ is a 
% projection onto $\text{Ran}(C^T)$, and $\mathbb{R}^m=\text{Ran}(C^T)\oplus\text{Ker}{(C)}$, we have 
% $E_0\in\text{Ker}(C)$. Since $E_0\in\text{Ker}(C)\cap\text{Ker}(\nabla^T)$
% we have from \eqref{identities} that
% %
% \begin{align}
% 0=C^TCE_0 
% = (-\nabla\nabla^T + I_d\otimes(\nabla^T\nabla))E_0
% = I_d\otimes(\nabla^T\nabla))E_0.
% \end{align}
% %
% From $\text{Ran}(\nabla)\subseteq\text{Ker}(C)$ in equation
% \eqref{subspaces} and $E_f=\nabla\varphi$ we have 
% $CE_f=C\nabla\varphi=0$. We also have from $\nabla^T\vecJ=0$ that
% $\vecJ\cdot E_f=\vecJ\cdot\nabla\varphi=\nabla^T\vecJ\cdot\varphi=0$.
% Finally, the volume average of $\nabla\varphi$ is a telescoping sum, 
% so $\langle E_f\rangle=0$. This concludes our proof of the theorem.
% \end{proof}


% \begin{corollary}
% %
% Let the electric and current fields $\vecE$ and $\vecJ$ satisfy
% %
% \begin{align}\label{eq:Maxwells_Equations_discrete}
% C\vecE=0, 
% \quad -\nabla^T J=0, 
% \quad J=\bsig E.
% \end{align}
% %
% Then, there exists a "potential" $\psi$ and vector $J_0$ such that
% $J=J_f+J_0$, where $J_f=C^T\psi$,
% $I_d\otimes(\nabla^T\nabla))J_0=0$, $\nabla^T J_f=0$, 
% $\langle\vecJ_f\cdot\vecE\rangle=0$, and $\langle\vecJ_f\rangle=0$.
% \end{corollary}
% %
% \begin{proof}
% From the resolution of the identity in equation \eqref{eq:res_of_identity}
% we have $J=(\Gamma_\times+\Gamma_0+\Gamma_\bullet)J$.
% Since $\Gamma_\times$ projects onto the $\text{Ran}(\nabla)$,
% $\mathbb{R}^m=\text{Ran}(\nabla)\oplus\text{Ker}(\nabla^T)$, and 
% $\vecJ\in\text{Ker}{\nabla^T}$ we have $\Gamma_\times J=0$. Denoting
% by $J_f=\Gamma_\bullet J$, since 
% $U_1^\bullet=C^T V_1^\bullet[\Sigma_1^\bullet]^{-1}$, we can write 
% $J_f=C^T\psi$, where
% $\psi=V_1^\bullet[\Sigma_1^\bullet]^{-1}[U_1^\bullet]^TJ$. Denote 
% $J_0=\Gamma_0 J$. Exactly the same as in the proof of Theorem
% \ref{thm:E_J_decomp}, we have $J_0= I_d\otimes(\nabla^T\nabla))J_0$. 
% From $\text{Ran}(C^T)\subseteq\text{Ker}(\nabla^T)$ in equation
% \eqref{subspaces} and $J_f=C^T\psi$ we have 
% $\nabla^T\vecJ_f=\nabla^TC^T\psi=0$. Exactly the same as in the proof 
% of Theorem \ref{thm:E_J_decomp}, we also have $\vecJ_f\cdot\vecE=0$ and 
% $\langle J_f\rangle=0$. This concludes our proof of the theorem.
% \end{proof}
% %
% We end this section by noting that in the full rank setting, where 
% $\Sigma$ has all strictly positive singular values, so $\Sigma_1=\Sigma$, 
% then 
% %
% \begin{align}\label{eq:Gamma}
% \Gamma_\times = \nabla(\nabla^T\nabla)^{-1}\nabla^T\,,
% \qquad
% \Gamma_\bullet = C^T(CC^T)^{-1}C\,.
% \notag
% \end{align}
% %
% The original formulations of this mathematical framework was given in terms of these projection matrices, or continuum generalizations \citep{Golden:CMP-473,Murphy:2015:CMS:13:4:825}. The formulation given in this section
% generalizes the discrete setting to cases where the matrix gradient is rank deficient,
% such as the case of periodic boundary conditions. This is necessary 


% \section{Proof of Theorem \ref{thm:Discrete_Spectral_Theorem_ACM_polycrystal}}
% \label{proof_thm}
% %
% In this section we provide the proof for Theorem
% \ref{thm:Discrete_Spectral_Theorem_ACM_polycrystal} and a projection
% method for a numerically efficient projection method for computation
% of spectral measures and effective parameters for uniaxial polycrystalline media.
% We will use the results from Section \ref{sec:field_decomp} but for notational 
% simplicity we will use $\Gamma$ instead of $\Gamma_\times$. Corollary
% \ref{cor:Discrete_Spectral_Theorem_ACM} below follows immediately from the proof of 
% Theorem \ref{thm:Discrete_Spectral_Theorem_ACM_polycrystal} and the results of Section
% \ref{sec:field_decomp}, which provides an integral representation
% for the effective resistivity $\brho^*$ involving the matrix $X_2\Gamma_\bullet X_2$ and is
% analogous to the representation of $\brho^*$ for the two-component composite 
% setting in \citep{Murphy:2015:CMS:13:4:825}. 

% From the close
% analogues of this polycrystalline setting with the two-component setting discussed in 
% Sec. \ref{sec:two_component}, the proof of this theorem is analogous to Theorem 2.1 in 
% \citep{Murphy:2015:CMS:13:4:825}. To shorten the theorem proof here, we will refer to 
% \citep{Murphy:2015:CMS:13:4:825} for some of the technical details.
% From Section \ref{sec:field_decomp} and the paragraph in \citep{Murphy:2015:CMS:13:4:825}
% containing equations 2.39 and 2.40,
% with $\chi_1$ replaced by $X_1$, we just need to prove that the functional
% $F_{jk}(s)=\langle(sI-X_1\Gamma X_1)^{-1}X_1\hat{e}_j\cdot\hat{e}_k\rangle$ 
% has the integral representation displayed in
% equation~\eqref{eq:Stieltjes_F_Discrete_polycrystal}. In the process, we will also establish a projection method for the numerically efficient, rigorous computation of
% $\mu_{jk}$. This projection method is summarized by
% equations~\eqref{Pi_coordinates_E}--\eqref{eq:Fs_W1} below.      





% In equation \eqref{eq:Projection_Matrices} we defined the real-symmetric 
% mutually orthogonal projection matrices $X_i$, $i=1,2$, in terms of the
% \emph{spatially dependent} rotation matrix $R$ and $C=\text{diag}(1,0,\ldots,0)$, 
% all matrices of size $d\times d$. The paragraph in \citep{Murphy:2015:CMS:13:4:825} 
% containing equations 2.28--2.30 describes how to bijectively map these $d\times d$
% matrices to $N\times N$ matrices that are \emph{not} spatially dependent, where $N=L^d d$.
% Under this mapping, $R$ becomes a banded rotation matrix satisfying
% $R^T=R^{-1}$ and $C$ becomes $C=\text{diag}(I_1,0_1,\ldots,0_1)$, where
% $I_1$ and $0_1$ are the identity and null matrices of size
% $N_1=L^d$, and the vector $e_1$ is mapped to $\hat{e}_1=(1,1,\ldots,1,0,0,\ldots,0)$,
% with $L^d$ ones in the first components and zeros in the rest of the components. 
% Writing $X_1\Gamma X_1=R^T(CR\,\Gamma \,R^TC)R$ we have 
% %
% \begin{align}\label{eq:Spec_Decomp_chi_Gamma_chi_Proof}
%   X_1\Gamma X_1&=
% R^T
% \left[
%   \begin{array}{ccc}
%     \Gamma_1&0_{12}\\
%     0_{21}&0_{22}   
%     \end{array}
% \right]
% R
% %\\
% %&=
% =
% R^T
% \left[
%   \begin{array}{ccc}
%     W_1\Lambda_1W_1^{\;T}&0_{12}\\
%     0_{21}&0_{22}
%     \end{array}
% \right]
% R
% \notag\\
% &=
% %=
% R^T
% \left[
%   \begin{array}{ccc}
%     W_1&0_{12}\\
%     0_{21}&I_2 
%     \end{array}
% \right]    
% \left[
%   \begin{array}{ccc}
%     \Lambda_1&0_{12}\\
%     0_{21}&0_{22}
%     \end{array}
% \right]    
% \left[
%   \begin{array}{ccc}
%     W_1^{\;T}&0_{12}\\
%     0_{21}&I_2
%     \end{array}
% \right]    
% R,
% %\notag
% \end{align}
% %
% where $I_2$ is the identity matrix of size $N_2\times N_2$, with
% $N_2=N-N_1=L^d(d-1)$, and $0_{ab}$ denotes a matrix of zeros of size
% $N_a\times N_b$, $a,b=1,2$. Moreover, $\Gamma_1$ is the $N_1\times N_1$ upper left principal
% sub-matrix of $R\,\Gamma\,R^T$ and $\Gamma_1=W_1\Lambda_1W_1^{\;T}$ is its eigenvalue
% decomposition. As $\Gamma_1$ 
% is a real-symmetric matrix, $W_1$ is an orthogonal
% matrix~\citep{Horn_Johnson-1990}. Also, since $R\,\Gamma\,R^T$ is a
% similarity transformation of a projection matrix and $C$ is a
% projection matrix, $\Lambda_1$ is a diagonal matrix with entries
% $\lambda_i^1\in[0,1]$, $i=1,\ldots,N_1$, along its
% diagonal~\citep{Horn_Johnson-1990,Demmel:1997}.     




% Consequently, equation~\eqref{eq:Spec_Decomp_chi_Gamma_chi_Proof}
% implies that the 
% eigenvalue decomposition of the matrix $X_1\Gamma X_1$ is given by 
% %$M=W\Lambda W^{\;T}$, where
% %
% \begin{align}\label{eq:Spec_Decomp_chi_Gamma_chi}
% X_1\Gamma X_1=W\Lambda W^{\;T},
% \qquad
% W=R^T\left[
%   \begin{array}{ccc}
%     W_1&0_{12}\\
%     0_{21}&I_2   
%     \end{array}
% \right],
% \quad
% \Lambda=\left[
%   \begin{array}{ccc}
%     \Lambda_1&0_{12}\\
%     0_{21}&0_{22}   
%     \end{array}
% \right].
% \end{align}
% %
% Here, $W$ is an orthogonal matrix satisfying $W^TW=WW^T=I$, $I$ is the
% identity matrix on $\mathbb{R}^N$, and $\Lambda$ is a diagonal matrix with
% entries $\lambda_i\in[0,1]$, $i=1,\ldots,N$, along its diagonal.




% The eigenvalue decomposition of the matrix $X_1\Gamma X_1$ in
% equation~\eqref{eq:Spec_Decomp_chi_Gamma_chi} demonstrates that its
% resolvent $(sI-X_1\Gamma X_1)^{-1}$ is well defined for all
% $s\in\mathbb{C}\backslash[0,1]$. In particular, by the orthogonality of the
% matrix $W$, it has the following useful representation
% $(sI-X_1\Gamma X_1)^{-1}=W(sI-\Lambda)^{-1}W^T$, where $(sI-\Lambda)^{-1}$ is a diagonal
% matrix with entries $1/(s-\lambda_i)$ along its diagonal. This, in
% turn, implies that the functional
% $F_{jk}(s)=\langle(sI-X_1\Gamma X_1)^{-1}X_1\hat{e}_j\cdot\hat{e}_k\rangle$ can be written as     
% %
% \begin{align}\label{eq:Matrix_Functional_proof}
%   F_{jk}(s)%=\langle(sI-X_1\Gamma X_1)^{-1}X_1\hat{e}_j\cdot\hat{e}_k\rangle
%           =\langle(sI-\Lambda)^{-1}\;[X_1W]^T\hat{e}_j\,\cdot\,W^T\hat{e}_k\rangle. 
% \end{align}
% %
% Since $R^T=R^{-1}$ and $X_1=R^TCR$,
% equation~\eqref{eq:Spec_Decomp_chi_Gamma_chi} implies that  
% %
% \begin{align}\label{eq:Projection_Eigenspace}
%   X_1W=R^T\left[
%   \begin{array}{ccc}
%     W_1&0_{12}\\
%     0_{21}&0_{22}  
%     \end{array}
% \right]
% \quad
% \Longrightarrow
% \quad
% X_1\vecw_i=
%   \begin{cases}
%   \vecw_i, &\text{~for~} i=1,\ldots,N_1,  \\
%   0,  &\text{~otherwise}.
% \end{cases}
% \end{align}
% %
% This and the formula for $W$ in~\eqref{eq:Spec_Decomp_chi_Gamma_chi}
% imply that 
% %
% \begin{align}\label{eq:Weights_chi}
%   [X_1W]^T\hat{e}_j\,\cdot\,W^T\hat{e}_k=[X_1W]^T\hat{e}_j\,\cdot\,[X_1W]^T\hat{e}_k.
% \end{align}
% %


% We are ready to provide the integral representation displayed
% in~\eqref{eq:Stieltjes_F_Discrete_polycrystal} for the functional $F_{jk}(s)$ in 
% equation~\eqref{eq:Matrix_Functional_proof}. Denote by
% $Q_i=\vecw_i\,\vecw_i^{\;T}$, $i=1,\ldots,N$, the symmetric, mutually
% orthogonal projection matrices, $Q_\ell\,Q_m=Q_\ell\,\delta_{\ell m}$, onto the eigen-spaces
% spanned by the orthonormal eigenvectors
% $\vecw_i$. Equation~\eqref{eq:Projection_Eigenspace} implies that 
% $X_1Q_i=Q_iX_1=X_1Q_iX_1$, as $X_1Q_i=Q_i$ for $i=1,\ldots,N_1$, $X_1Q_i=0$
% otherwise, and $X_1$ is a symmetric matrix. These properties allow us to write
% the quadratic form $[X_1W]^T\hat{e}_j\,\cdot\,[X_1W]^T\hat{e}_k$ as  
% %
% \begin{align}\label{eq:Quadratic_form}
%   [X_1W]^T\hat{e}_j\,\cdot\,[X_1W]^T\hat{e}_k =
%   W^T\hat{e}_j\,\cdot\,W^T\hat{e}_k =
%   \sum_{i=1}^N
%   (\vecw_i\cdot\hat{e}_j)(\vecw_i\cdot\hat{e}_k) =
%   \sum_{i=1}^N Q_i\hat{e}_j\cdot\hat{e}_k =
%  \sum_{i=1}^N X_1Q_i\hat{e}_j\cdot\hat{e}_k\,.
% \end{align}
% %
% This and equations~\eqref{eq:Matrix_Functional_proof}
% and~\eqref{eq:Weights_chi} yield 
% %
% \begin{align}\label{eq:Stieltjes_F_DiscretE}
%   F_{jk}(s)=\int_0^1\frac{\d\mu_{jk}(\lambda)}{s-\lambda}\,, \quad
%   \d\mu_{jk}(\lambda)=\sum_{i=1}^N\langle \delta_{\lambda_i}(\d\lambda)\;X_1\,Q_i\hat{e}_j\cdot\hat{e}_k\rangle\,.
% \end{align}
% %
% This 
% concludes our proof of Theorem~\ref{thm:Discrete_Spectral_Theorem_ACM_polycrystal}


% \begin{corollary}

% \label{cor:Discrete_Spectral_Theorem_ACM}
% For each $\omega\in\Omega$, let $M(\omega)=W(\omega)\Lambda(\omega)\,W(\omega)$ be 
% the eigenvalue decomposition of the real-symmetric matrix $M(\omega)=X_2(\omega)\,\Gamma_\bullet\,X_2(\omega)$. Here, the columns of the matrix $W(\omega)$
% consist of the orthonormal eigenvectors $\vecwt_i(\omega)$, $i=1,\ldots,N$,
% of $M(\omega)$ and the diagonal matrix $\Lambda(\omega)={\rm diag}(\lambda_1(\omega),\ldots,\lambda_N(\omega))$
% involves its eigenvalues $\lambda_i(\omega)$. 
% Denote $Q_i=\vecwt_i\,\vecwt_i^{\;T}$ the projection matrix onto the eigen-space
% spanned by $\vecwt_i$.
% The curent field $\vecJ(\omega)$ satisfies $\vecJ(\omega)=\vecJ_0+\vecJ_f(\omega)$, with
% $\vecJ_0=\langle\vecJ(\omega)\rangle$ and $\Gamma_\bullet\vecJ(\omega)=\vecJ_f(\omega)$, 
% and the effective complex resistivity tensor $\brho^*$ has components
% $\rho^*_{jk}$,
% %   $\sigma_{jk}^*=\sigma_2\,m_{jk}(h)$, 
% $j,k=1,\ldots,d$,  which satisfy       
% %
% \begin{align}\label{eq:Stieltjes_F_Discrete_rho}
%   &\rho^*_{jk}=\sigma_1^{-1}(\delta_{jk}-E_{jk}(s)), 
%   &&E_{jk}(s)=\int_0^1\frac{\d\eta_{jk}(\lambda)}{s-\lambda}\,, 
%   &&\d\eta_{jk}(\lambda)=\sum_{i=1}^N\langle \delta_{\lambda_i}(\d\lambda)\;X_2\,Q_i\hat{e}_j\cdot\hat{e}_k\rangle\,.  
% \end{align}
% % 
% \end{corollary}


% \section{Projection method}\label{projection_method}


% In this section we provide a formulation for a numerically efficient 
% \emph{projection method} for computation of spectral measures and 
% effective parameters for uniaxial polycrystalline media. Note that the 
% sum inequation~\eqref{eq:Stieltjes_F_DiscretE} runs only over the index 
% set $i=1,\ldots,N_1$, as equation~\eqref{eq:Projection_Eigenspace} implies
% that the masses $X_1\,Q_i\hat{e}_j\cdot\hat{e}_k$ of the measure $\mu_{jk}$
% are zero for $i=N_1+1,\ldots,N$.  Denote by $\lambda_i^1$ and $\vecw_i^{\,1}$,
% $i=1,\ldots,N_1$, the eigenvalues and eigenvectors of the $N_1\times N_1$
% matrix $\Gamma_1=W_1\Lambda_1W_1^{\;T} $, defined in
% equation~\eqref{eq:Spec_Decomp_chi_Gamma_chi_Proof}. Now, write        
% %
% \begin{align}\label{Pi_coordinates_E} 
% R\hat{e}_j=
%   \left[
%   \begin{array}{ccc}
%     \hat{e}_j^{\,r_1}\\
%     \hat{e}_j^{\,r_2}
%     \end{array}
% \right],
% \end{align}
% %
% where $\hat{e}_j^{\,r_1}\in\mathbb{R}^{N_1}$ and
% $\hat{e}_j^{\,r_2}\in\mathbb{R}^{N_2}$. Therefore, writing the
% matrix $X_1W$ in equation~\eqref{eq:Projection_Eigenspace} in block
% diagonal form, $X_1W=R^T\text{diag}(W_1,0_{22})$, we have that
% %
% \begin{align}\label{eq:Reduced_Weights}
%   [X_1W]^T\hat{e}_j\,\cdot\,[X_1W]^T\hat{e}_k=[\text{diag}(W_1^T,0_{22})R\hat{e}_j]
%                                     \cdot[\text{diag}(W_1^T,0_{22})R\hat{e}_k]
%                                    =[W_1^T\hat{e}_j^{\,r_1}]
%                                    \cdot[W_1^T\hat{e}_k^{\,r_1}]. 
% \end{align}
% %
% Denote by $Q^1_i=\vecw^{\,1}_i[\vecw^{\,1}_i]^{\;T}$, $i=1,\ldots,N_1$,
% the mutually orthogonal projection matrices,
% $Q^1_\ell\,Q^1_m=Q^1_\ell\,\delta_{\ell m}$, onto the eigen-spaces spanned by the
% orthonormal eigenvectors
% $\vecw^{\,1}_i$. Equations~\eqref{eq:Matrix_Functional_proof},~\eqref{eq:Weights_chi},   
% and~\eqref{eq:Reduced_Weights} then yield    
% %
% \begin{align}\label{eq:Fs_W1}
%   F_{jk}(s)=\langle(sI_1-\Lambda_1)^{-1}[W_1^T\hat{e}_j^{\,r_1}]
%                        \cdot[W_1^T\hat{e}_k^{\,r_1}]\rangle
%           =\left\langle\sum_{i=1}^{N_1} 
%           \frac{Q^1_i\hat{e}_j^{\,r_1}\cdot\hat{e}_k^{\,r_1}                
%               }{s-\lambda_i^1}
%               \right\rangle.        
% \end{align}
% %



% Equation~\eqref{eq:Fs_W1} demonstrates that only the spectral
% information of the matrices $W_1$ and $\Lambda_1$ contribute to the functional
% representation for $F_{jk}(s)$ in~\eqref{eq:Matrix_Functional_proof}
% and its integral representation
% in~\eqref{eq:Stieltjes_F_Discrete_polycrystal}. From a computational standpoint,  
% this means that only the eigenvalues and eigenvectors of the $N_1\times N_1$
% matrix $\Gamma_1$ need to be computed in order to compute the spectral
% measures underlying the integral representations of the effective
% parameters for finite lattice systems. This is extremely cost
% effective as the numerical cost of finding all
% the eigenvalues and eigenvectors of a real-symmetric $N\times N$ matrix is
% $O(N^3)$~\citep{Demmel:1997}, so $N_1=N/d$ implies the computational cost
% of the projection method is reduced by a factor of $d^3$. 


\noappendix       %% use this to mark the end of the appendix section. Otherwise the figures might be numbered incorrectly (e.g. 10 instead of 1).

%% Regarding figures and tables in appendices, the following two options are possible depending on your general handling of figures and tables in the manuscript environment:

%% Option 1: If you sorted all figures and tables into the sections of the text, please also sort the appendix figures and appendix tables into the respective appendix sections.
%% They will be correctly named automatically.

%% Option 2: If you put all figures after the reference list, please insert appendix tables and figures after the normal tables and figures.
%% To rename them correctly to A1, A2, etc., please add the following commands in front of them:

% \appendixfigures  %% needs to be added in front of appendix figures

% \appendixtables   %% needs to be added in front of appendix tables

%% Please add \clearpage between each table and/or figure. Further guidelines on figures and tables can be found below.


\authorcontribution{All authors contributed to planning the work and writing and editing the manuscript.} %% this section is mandatory
 

\competinginterests{The authors declare that they have no conflict of interest.} %% this section is mandatory

%\disclaimer{TEXT} %% optional section

\begin{acknowledgements}
We gratefully acknowledge support from the 
Applied and Computational Analysis Program and
the Arctic and Global Prediction Program 
at the US Office of Naval Research
through grants N00014-18-1-2552, ONR Grant N00014-21-1-2909,
N00014-13-10291, 
N00014-15-1-2455, and N00014-18-1-2041.
We are also grateful 
for support from the
Division of Mathematical 
Sciences and the Division of Polar Programs at the US National 
Science Foundation (NSF) through Grants DMS-0940249,  %mcrn
DMS-2136198, %RTG
DMS-2111117,
DMS-2206171,
DMS-1715680, %rmt
and DMS-1413454.
Finally, we would like to
thank the NSF Math Climate Research Network (MCRN),
and especially Chris Jones,
for supporting this work.
\end{acknowledgements}




%% REFERENCES

%% The reference list is compiled as follows:

%%% nbenjamin murphy Using bib files
% % \small
\clearpage
% \pagestyle{empty}
%\setcounter{page}{1}

%\bibliographystyle{unsrt}
%\bibliographystyle{unsrtnat}
% \bibliographystyle{plain}
%\bibliographystyle{./ametsoc}
\vskip -15cm
%\bibliography{jabref_database,strong_results,golden,cmg,hyde,opp08,Murphy_DMS_2013,elastic}

%\bibliography{golden,newgolden,murphy,elastic,WaveRef}

\bibliographystyle{copernicus}
\bibliography{golden.bib,newgolden.bib,murphy.bib,elastic.bib,WaveRef.bib}
%\begin{thebibliography}{}
%
%\bibitem[AUTHOR(YEAR)]{LABEL1}
%REFERENCE 1
%
%\bibitem[AUTHOR(YEAR)]{LABEL2}
%REFERENCE 2
%
%\end{thebibliography}

%% Since the Copernicus LaTeX package includes the BibTeX style file copernicus.bst,
%% authors experienced with BibTeX only have to include the following two lines:
%%
%% \bibliographystyle{copernicus}
%% \bibliography{example.bib}
%%
%% URLs and DOIs can be entered in your BibTeX file as:
%%
%% URL = {http://www.xyz.org/~jones/idx_g.htm}
%% DOI = {10.5194/xyz}


%% LITERATURE CITATIONS
%%
%% command                        & example result
%% \citet{jones90}|               & Jones et al. (1990)
%% \citep{jones90}|               & (Jones et al., 1990)
%% \citep{jones90,jones93}|       & (Jones et al., 1990, 1993)
%% \citep[p.~32]{jones90}|        & (Jones et al., 1990, p.~32)
%% \citep[e.g.,][]{jones90}|      & (e.g., Jones et al., 1990)
%% \citep[e.g.,][p.~32]{jones90}| & (e.g., Jones et al., 1990, p.~32)
%% \citeauthor{jones90}|          & Jones et al.
%% \citeyear{jones90}|            & 1990



%% FIGURES

%% When figures and tables are placed at the end of the MS (article in one-column style), please add \clearpage
%% between bibliography and first table and/or figure as well as between each table and/or figure.

% The figure files should be labelled correctly with Arabic numerals (e.g. fig01.jpg, fig02.png).


%% ONE-COLUMN FIGURES

%%f
%\begin{figure}[t]
%\includegraphics[width=8.3cm]{FILE NAME}
%\caption{TEXT}
%\end{figure}
%
%%% TWO-COLUMN FIGURES
%
%%f
%\begin{figure*}[t]
%\includegraphics[width=12cm]{FILE NAME}
%\caption{TEXT}
%\end{figure*}
%
%
%%% TABLES
%%%
%%% The different columns must be seperated with a & command and should
%%% end with \\ to identify the column brake.
%
%%% ONE-COLUMN TABLE
%
%%t
%\begin{table}[t]
%\caption{TEXT}
%\begin{tabular}{column = lcr}
%\tophline
%
%\middlehline
%
%\bottomhline
%\end{tabular}
%\belowtable{} % Table Footnotes
%\end{table}
%
%%% TWO-COLUMN TABLE
%
%%t
%\begin{table*}[t]
%\caption{TEXT}
%\begin{tabular}{column = lcr}
%\tophline
%
%\middlehline
%
%\bottomhline
%\end{tabular}
%\belowtable{} % Table Footnotes
%\end{table*}
%
%%% LANDSCAPE TABLE
%
%%t
%\begin{sidewaystable*}[t]
%\caption{TEXT}
%\begin{tabular}{column = lcr}
%\tophline
%
%\middlehline
%
%\bottomhline
%\end{tabular}
%\belowtable{} % Table Footnotes
%\end{sidewaystable*}
%
%
%%% MATHEMATICAL EXPRESSIONS
%
%%% All papers typeset by Copernicus Publications follow the math typesetting regulations
%%% given by the IUPAC Green Book (IUPAC: Quantities, Units and Symbols in Physical Chemistry,
%%% 2nd Edn., Blackwell Science, available at: http://old.iupac.org/publications/books/gbook/green_book_2ed.pdf, 1993).
%%%
%%% Physical quantities/variables are typeset in italic font (t for time, T for Temperature)
%%% Indices which are not defined are typeset in italic font (x, y, z, a, b, c)
%%% Items/objects which are defined are typeset in roman font (Car A, Car B)
%%% Descriptions/specifications which are defined by itself are typeset in roman font (abs, rel, ref, tot, net, ice)
%%% Abbreviations from 2 letters are typeset in roman font (RH, LAI)
%%% Vectors are identified in bold italic font using \vec{x}
%%% Matrices are identified in bold roman font
%%% Multiplication signs are typeset using the LaTeX commands \times (for vector products, grids, and exponential notations) or \cdot
%%% The character * should not be applied as mutliplication sign
%
%
%%% EQUATIONS
%
%%% Single-row equation
%
%\begin{equation}
%
%\end{equation}
%
%%% Multiline equation
%
%\begin{align}
%& 3 + 5 = 8\\
%& 3 + 5 = 8\\
%& 3 + 5 = 8
%\end{align}
%
%
%%% MATRICES
%
%\begin{matrix}
%x & y & z\\
%x & y & z\\
%x & y & z\\
%\end{matrix}
%
%
%%% ALGORITHM
%
%\begin{algorithm}
%\caption{...}
%\label{a1}
%\begin{algorithmic}
%...
%\end{algorithmic}
%\end{algorithm}
%
%
%%% CHEMICAL FORMULAS AND REACTIONS
%
%%% For formulas embedded in the text, please use \chem{}
%
%%% The reaction environment creates labels including the letter R, i.e. (R1), (R2), etc.
%
%\begin{reaction}
%%% \rightarrow should be used for normal (one-way) chemical reactions
%%% \rightleftharpoons should be used for equilibria
%%% \leftrightarrow should be used for resonance structures
%\end{reaction}
%
%
%%% PHYSICAL UNITS
%%%
%%% Please use \unit{} and apply the exponential notation


\end{document}
