\documentclass[12pt]{article}





\usepackage{graphicx}
\usepackage{amssymb}   
\usepackage{amsmath}
\usepackage{amsfonts}
%\usepackage{dcolumn}



\usepackage{natbib}
\bibpunct{(}{)}{,}{a}{}{;}






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

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








%\setlength{\topmargin}{-0.5in}
%\setlength{\headheight}{0in}
%\setlength{\headsep}{0in}
%\oddsidemargin 0in
%\textheight 9in
%\textwidth 6.5in





%\setlength{\topmargin}{-1.0in}
%\setlength{\headheight}{0in}
%\setlength{\headsep}{0in}
%\oddsidemargin -0.35in
%\textheight 8in
%\textwidth 6in
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\half}{\frac{1}{2}}
\newcommand{\tf}{\tilde{f}}






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








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











\newcommand{\vv}{\textbf{v}}
\newcommand{\xx}{\textbf{x}}
\newcommand{\yy}{\textbf{y}}
\newcommand{\uu}{\textbf{u}}
\newcommand{\UU}{\textbf{U}}
\newcommand{\ff}{\textbf{f}}
\newcommand{\KK}{\textbf{K}}
\newcommand{\VV}{\textbf{V}}
\newcommand{\kk}{\textbf{k}}
\newcommand{\ww}{\textbf{w}}
\newcommand{\II}{\textbf{I}}

\newcommand\bsig{\mbox{\boldmath${\sigma}$}}
\newcommand\beps{\mbox{\boldmath${\epsilon}$}}
\newcommand\bnab{\mbox{\boldmath${\nabla}$}}
\newcommand\bcdot{\mbox{\boldmath${\cdot}$}}
\newcommand\btimes{\mbox{\boldmath${\times}$}}
\newcommand\bxi{\mbox{\boldmath${\xi}$}}  







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







\newcommand{\hf}{\frac{1}{2}}





\begin{document}




\begin{center}


{\LARGE {\bf Spectral Analysis and Connectivity of  \\
   \vspace{1ex}
    Porous Microstructures in Bone} }


%\author


\vspace{10ex}



{\large Kenneth M. Golden, N. Benjamin Murphy, and Elena Cherkaev


\vspace{8ex}


University of Utah, Department of Mathematics \\
\vspace{4ex}
155 S 1400 E RM 233, Salt Lake City, UT, 84112-0090 USA}




\vspace{20ex}

Corresponding Author: Kenneth M. Golden\\

\vspace{4ex}

Phone: 801-581-6176 $\;\;\;\;$
Fax: 801-581-4148   $\;\;\;\;$
Email: golden@math.utah.edu



\vspace{8ex}



\noindent
Keywords:
bone, sea ice, porosity, percolation, spectral reconstruction 



\vspace{8ex}



\noindent
Approximate word count: 3,800






\end{center}


\baselineskip32pt

%\baselineskip14pt

%\maketitle 

\eject

\begin{center}


{\bf Abstract}

\end{center} 


\vspace{3ex}

%\begin{abstract}
\noindent
Cancellous bone is a porous composite of calcified 
tissue interspersed with soft marrow.
Sea ice is also a porous composite, consisting of pure ice with 
brine, air, and salt inclusions.
%whose volume fraction and connectivity
%vary significantly with temperature.
Interestingly, the microstructures 
of bone and sea ice exhibit notable similarities.
In recent years, we have developed mathematical and
experimental techniques for imaging and characterizing the 
brine microstructure of sea ice, such as its volume
fraction and connectivity, as well as
a range of theoretical approaches for studying
fluid, thermal, and electromagnetic transport in sea ice.
Here we explore the application of our sea ice techniques 
to investigate trabecular bone. For example, 
percolation theory that
quantifies brine connectivity and its thermal evolution 
can also help assess the impact of osteoporosis on trabecular 
structure. Central to our approach is the 	
{\it spectral measure} of a composite material, which
contains detailed information about the mixture geometry,
and can be used in powerful integral representations
to compute the effective properties. The spectral measure 
is obtained from the eigenvalues and eigenvectors
of a self adjoint operator determined exclusively by the 
composite microgeometry. Here we compute the spectral measures
for discretizations of images of healthy and osteoporotic bone.
The measures are used to compute the effective 
electromagnetic properties of the bone specimens. 
These data are then inverted to reconstruct the porosity 
of the original specimens, with excellent agreement.
%\end{abstract}




\vspace{5ex}













\eject





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%









\section{Introduction}




Bone displays a complex, porous microstructure whose
characteristics depend on its macrostructure, whether
cortical or cancellous, as well as age and
health of the individual. 
The strength of bone and its ability to resist
fracture depend strongly on this porous microstructure,
and in particular, on the {\it quality} of the 
connectedness of the hard, solid phase
\citep{Odgaard:B:1997,Kabel:B:1999,Hollister:JB:1994}.
For example, in dense cortical bone the pores can be sparse and 
disconnected, yet exhibit increasing volume fraction and 
connectivity with the onset of osteoporosis.
Cancellous bone displays a broad range of 
biconnected microstructures,
ranging from a solid network of connected
trabeculae containing numerous 
connected pores, to sparse
solid fibers within a 
dominant, connected pore space. 
With the onset of osteoporosis,
cancellous bone can become
more disconnected and remaining connections can
become more tenuous or fragile. 
There have been many studies of bone structure
(including porosity) and mechanics, 
and how they depend on aging and other factors
\citep{Ural:JB:2006,Fritsch:JTB:2007,Nalla:B:2004,%
Hildebrand:JBMR:1999,Coelho:JB:2009,Burghardt:JBMR:2010,%
Hollister:JB:1994}.






In this paper we consider what
percolation theory, or the mathematical theory
of connectedness, can tell us about bone structure.
In particular, we investigate the spectral 
measures for two examples of 
porous bone microstructure. A spectral measure
is the key mathematical object appearing in
integral representations for effective transport
and elastic properties of two phase composites
\citep{Golden:CMP-473,Kantor:JMPS:1984}.
It contains, in principle, most of the geometrical 
information about the composite. For example,
the mass of the measure is the porosity, or
volume fraction of one phase, and
the absence of large-scale connectivity is
associated with a gap in the spectrum.
Here we intend
to introduce methods of analyzing bone microstructure
and its properties which may eventually help 
in clinical applications.
Steps in this direction have been taken already in 
\citep{Bonifasi:2006,Bonifasi:2008,Bonifasi:2009,Bonifasi:2009a},
where {\it inverse homogenization theory} 
\citep{Cherkaev:2001,Cherkaeva:WRM-437}
has been extended to the evaluation of bone structure,
which is discussed in detail below.






Here we observe that some of the natural questions
we ask about sea ice, as well as the methods to 
address them, are also relevant for studies
of bone. 
For sea ice 
we have applied percolation theory to understand
the thermal evolution of the brine microstructure
and its connectedness
\citep{Golden:S-2238,Golden:2007:GRL,Pringle:JGR:2009}.
We have also developed X-ray CT imaging methods
and pore structure analysis 
\citep{Golden:2007:GRL,Pringle:JGR:2009}
to investigate
the microstructure and verify our earlier conjecture
of a connectedness and fluid transport transition at a brine volume fraction
of about 5\%, which corresponds to a critical 
temperature of $-5^{\circ}$C for a typical bulk
salinty of 5 parts per thousand,
known as the {\it rule of fives}
\citep{Golden:S-2238}.
In studying the role of sea ice
in the climate system, we are then interested in how
this complex microstructure determines the 
fluid \citep{Golden:2007:GRL}
and thermal transport properties
of sea ice.
The electromagnetic properties of sea ice with applications
to remote sensing are studied extensively in 
\citep{Golden:JGR-13699,Golden:IMA-97,Golden:TGRS-1655,%
Golden:TGRS-1675,Gully:2007:Physica}. An overview
is given by \cite{Golden:NAMS:2009}. 
Microstructural imaging and characterization of bone 
\citep{Fritsch:JTB:2007,Hildebrand:JBMR:1999,Kazakia:JBMR:2008,%
Issever:JBMR:2009,Burghardt:JBMR:2010}
is fundamental to studies of its effective 
mechanical properties
\citep{Ural:JB:2006,Nalla:B:2004,Oyen:JB:2008,%
Coelho:JB:2009,Hollister:JB:1994},
as well as its fluid
\citep{Piekarski:N:1977,KnotheTate:JB:2008,KnotheTate:BM:2009},
electromagnetic \citep{Sierpowska:PMB:2003,Singh:CORR:1984,%
Saha:ITBE:1992,Williams:ABE:1996}, and thermal
\citep{Davidson:MEP:2000,Biyikli:JBMR:1986} transport properties.


It should be mentioned that there have been other mathematical
approaches to characterizing the connectivity of
trabecular architecture and its relation to bone quality
and strength. In particular, there has been substantial
research focused on using the so-called 
{\it Euler-Poincar\'{e} characteristic},
an integral geometrical
quantity, to provide an estimate of the connectivity of the
pore space structure
\citep{Odgaard:B:1993,Odgaard:B:1997,Kabel:B:1999,Jinnai:B:2002,Roque:RBR:2009}.
An important aspect of this characteristic is that
it does not change under deformation or scaling of an object --
it is a topological invariant which has been investigated
to distinguish between 
osteoporotic and healthy bone.
The percolation and spectral characterization of 
bone microstructure presented here is more statistical 
in nature, rooted in issues of long range
order and connectedness, and motivated by statistical
mechanics and the theory of random composites.
Understanding the relations between
our approach and the Euler-Poincar\'{e} characteristic
is a topic for future research.





\section{Mathematical Methods}


\subsection{Analytic continuation and the spectral measure}

{\it Homogenization} denotes a field of applied mathematics
where the goal is to find a homogeneous medium which behaves
macroscopically just like a given inhomogeneous medium.
We briefly describe the analytic continuation method
for studying the effective properties of composite materials
\citep{Bergman:PRL-1285,Milton:APL-300,Golden:CMP-473,Golden:IMA-97},
or for homogenizing inhomogeneous media.
Let $\epsilon(x,\o)$ be a spatially stationary random field
in $x\in \mathbb{R}^d$, where $d$ is spatial dimension,
and $\o \in \O$, where $\O$ is the
set of all realizations of the random medium,
which represents the local values 
of the complex permittivity. This key parameter is defined in 
basic electromagnetism and determines the propagation
properties of an electromagnetic wave in a medium.
The real part represents the polarizability of the medium,
and the imaginary part the losses. 
We assume  $\epsilon(x)$ (where we supress the $\o$ notation)
takes the values 
$\epsilon_1$ in bone marrow
and $\epsilon_2$ in bone, where both of these 
parameters depend on frequency,
and write
%\begin{equation}                                 \label{eq2.1}
$\epsilon(x) = \epsilon_1 (1-\chi (x)) + \epsilon_2 \chi(x),$
%\end{equation}
where $\chi$ is the characteristic function of the bone, 
which equals one for all realizations $\o \in \O$
having bone at $x$, and equals zero otherwise.
Let $E (x)$ and $D(x)$ be the stationary
random electric and displacement fields satisfying
the constitutive law
%\begin{equation}                                 \label{eq2.2}
$    D (x)= \epsilon (x) E (x) $
%\end{equation}
and the equations
\begin{equation}                                 \label{eq2.3}
\nabla \cdot D  = 0, \;\;\;\;\;
\nabla \times E  = 0,
\end{equation}
with
%\begin{equation}                                   \label{eq2.5}
$\langle E (x) \rangle = e_k,$
%\end{equation}
where $e_k$ is a unit vector in the $k^{th}$ direction for
some $k = 1, \ldots, d$, and $\langle \cdot \rangle$ means
an ensemble average over $\O$ or spatial average over all of
$\mathbb{R}^d$.

The effective complex permittivity tensor ${\beps}^*$
is defined by
\begin{equation}                                    \label{eq2.6}
    \langle D \rangle= {\beps}^* \langle E \rangle.
\end{equation}
Without loss of generality, 
we focus on one diagonal coefficient
$\e^*~=~\e^*_{kk}$. Due to homogeneity,
$\e^*(a \epsilon_1, a \epsilon_2)= a
\e^*(\epsilon_1,\epsilon_2)$, for any complex number $a$,
$\e^*$ depends only on the ratio
$h = \epsilon_2 / \epsilon_1$, and we define
$m(h) =\e^* / \epsilon_1$.
The two main properties of $m(h)$ are that it is analytic
off $(- \infty, 0]$ in the $h$--plane, and that it
maps the
upper half plane to the upper half plane,
%\citep{Bergman:PRC-377,Golden:CMP-473},
so that it is an example of a Herglotz, or Stieltjes function.
The key step in the
analytic continuation method is obtaining
an integral representation for $\e^*$.

It is more convenient to work with 
the function %\citep{Bergman:PRC-377}
%\begin{equation}                                    \label{eq2.7}
$F(s) = 1 - m(h)$, where
$s = 1/(1-h)$,
%\end{equation}
which is analytic off $[0,1]$ in the $s$--plane.  It
was proven \citep{Golden:CMP-473,Bergman:PRC-377} 
that $F(s)$ has the 
representation
\begin{equation}                                     	\label{eq2.8}
    F(s)=\int_0^1 {\frac{d\mu(\lambda)}{s-\lambda}}\;, 
\qquad  s\notin[0,1]\;,
\end{equation}
where $\mu$ is a positive measure on $[0,1]$. 
This formula arises from the resolvent
representation  of the electric field 
$E = s(s + \Gamma \chi)^{-1}e_k$,
where $\Gamma  = \nabla (- \Delta )^{-1} \nabla \cdot$ 
and $\Delta=\nabla^2$ is the Laplacian, 
yielding
\begin{equation}                              	\label{gold1.eq.3a}
F(s) = \langle \chi [(s + \Gamma \chi)^{-1} e_k] \cdot e_k \rangle.
\end{equation}
%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)$ with weight $\chi$ in
the inner product, $\Gamma \chi$ is a bounded self adjoint
operator. % with $\parallel \Gamma \chi_1 \parallel \leq 1$.
Formula (\ref{eq2.8})
is the spectral representation of the resolvent, 
and $\mu$ is a {\it spectral measure} of $\Gamma \chi$, in the 
$e_k$ state.
The integral representation
(\ref{eq2.8}) separates the
parameter information in $s$
from information about the mixture geometry 
contained in $\mu$.
Statistical assumptions about the 
geometry are incorporated into $\mu$ via its moments
%\begin{equation}                                   \label{eq2.9}
$\mu_n = \int^1_0 \lambda^n d \mu (\lambda)$,
which can be calculated from the
correlation functions of the random medium, with 
\begin{equation}
\mu_n = (-1)^n \langle \chi 
[( \Gamma \chi )^n e_k ] \cdot e_k \rangle.
\end{equation}
The mass $\mu_0$ is the 
bone volume fraction $\phi$ (with $1-\phi$ the porosity),
\begin{equation}                                   \label{BVF}
\mu_0 = \int^1_0  d \mu (z) = 
\langle \chi \rangle = \phi.
\end{equation}
%In the notation we use here, the porosity of the
%bone composite, or volume fraction of the marrow phase,
%is $1-\phi$.


A principal application of the analytic continuation method
is to derive {\it forward bounds} on $\e^*$ given partial
information on the microgeometry
\citep{Bergman:PRL-1285,Milton:APL-300,%
Golden:IMA-97,Golden:TGRS-1655,Gully:2007:Physica}.
The objective of {\it inverse bounds}
is to use data about 
the electromagnetic response of a composite to bound 
its structural parameters such as $\phi$ 
\citep{Cherkaeva:WRM-437,Golden:TGRS-1675,%
McPhedran:1982,McPhedran:1990,Cherkaev:1996,Tripp:1998}.
Below we will use a more sophisticated inversion
scheme over frequency
to reconstruct the spectral measure, and obtain
accurate estimates of bone volume
fraction $\phi$.
A powerful application of this inverse 
homogenization approach \citep{Cherkaev:2001}   
is that once the spectral measure has been 
reconstructed, it can be used to calculate other effective 
properties such as DC electrical conductivity or thermal 
conductivity, which have
%These other transport properties have 
analogous integral representations. 
%but involving the 
%ratio of the relevant component parameters. 
Thus, through this spectral coupling, other transport properties can 
be estimated indirectly from EM monitoring 
\citep{Cherkaev:2001,Cherkaev:2003,Cherkaev_Zhang:2003}.

 





















 




\subsection{Reconstruction of the spectral measure}



The problem of extracting microstructural information
from bulk property measurments %using analytic continuation
was introduced previously by 
\cite{McPhedran:1982} and \cite{McPhedran:1990}.  
They estimated the
volume fraction of one component in a two phase mixture from
measurements of the effective complex permittivity of the composite
material. An analytical approach
to estimating the volume fractions of materials in a composite was
developed \citep{Cherkaev:1996,Tripp:1998},
and extended to the problem of finding bounds on the
microstructural parameters for isotropic composite materials
in \citep{Cherkaeva:WRM-437}. 
The analytical method is based on the integral
representation in (\ref{eq2.8}) and gives explicit formulas 
for bounds on the volume fraction of one 
of the constituents \citep{Cherkaeva:WRM-437}.












The problem of characterizing
microstructural information from effective property measurements was
formulated as an inverse problem for the spectral measure $\mu$ in
the Stieltjes analytic representation \citep{Cherkaev:2001}. 
Uniqueness of the reconstruction
of the spectral measure \citep{Cherkaev:2001} gives a
basis for the theory of inverse homogenization. 
The spectral measure $\mu$ in the
integral representation (or spectral function) 
can be uniquely reconstructed if
measurements of the effective properties of the composite are
available along an arc in the complex $s$--plane 
\cite{Cherkaev:2001}.
Such data can be
obtained from measurements in an interval of frequency provided that
at least one of the constituents is frequency dependent.









The numerical problem of reconstructing the spectral function is
ill-posed, that is, the solution is highly sensitive to changes in 
the data,
and requires regularization \citep{Cherkaev:2001}. Several
regularized algorithms have been developed
\citep{Cherkaev:2001,Cherkaev:2003,%
Cherkaev_Zhang:2003}. To obtain a stable
reconstruction of the spectral measure, an inversion method based on
constrained rational approximation of the spectral function was
developed in \citep{Zhang:2008,Zhang:2009}.




The inverse homogenization method was extended to the viscoelastic problem in
\citep{Bonifasi:2006,Bonifasi:2008} for torsion
of a cylinder whose microstructure does not change in the axial direction. With this
simplified model of bone, this method was used successfully to recover porosity
from measurements of the effective shear modulus simulated using micro-CT images
of cancellous bone \citep{Bonifasi:2008, Bonifasi:2009a}.
An electrical impedance spectroscopy method developed in
\citep{Bonifasi:2009} shows potential for reconstructing bone
porosity from effective electrical measurements.
The method is based on reconstruction of the spectral 
function from electrical measurements and calculating the
zeroth moment, which gives the volume fraction of one of 
the bone components \citep{Bonifasi:2009}. 
The method was applied to the
effective electrical properties of cancellous bone, and
numerically calculated using micro-CT images of human 
vertebrae to find the bone porosity. 
It was shown that bone porosity can be accurately 
reconstructed even in the presence of errors in the 
measured data and uncertainty in the properties of bone 
and bone marrow.











 The problem of reconstructing the
spectral measure $\mu$  can be reduced to an inverse potential
problem. It was  shown in \citep{Cherkaev:2001} that 
$F(s)$ admits a representation
as  a logarithmic potential of the measure $\mu$:
\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 potential $\Phi$ solves the Poisson equation
$ - \Delta \Phi = \rho$,  
where  $\rho(\lambda)$  is  a density on
$[0,1]$. A solution to the forward problem is given by the
Newtonian potential with $ \mu(d\lambda) = \rho(\lambda)d\lambda$. 
The inverse problem is to find $\rho(\lambda)$ (or $\mu$) given
values of $\Phi, \partial \Phi \, /
\partial n $,  or  $\nabla \Phi$.



%It is known that the reconstruction problem for the logarithmic
%potential is extremely ill-posed  and  requires  regularization to
%develop a stable numerical algorithm. 
To construct the solution to the inverse problem we
formulate the related minimization problem:
\begin{equation}                               \label{minimization}
\mbox{min}_{\mu} \; || A \mu - F || \,,
\end{equation}
where $A\mu = 
\frac{\D}{\D s}\int_0^1 \ln\, | s - \lambda | \, \,d \mu(\lambda) $, 
%is the integral operator in  (\ref{ln}), 
the norm  is
the $L^2-$norm, and $F$ is the given function of the measured data, $F(s)
= 1-\e^*(s)/\e_1, ~ s \in {\mathbb{C}} $.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The solution of the 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 a
regularization technique.

A regularization algorithm developed in \citep{Cherkaev:2001} is
based on constrained minimization.  It  introduces   a stabilization
functional   which  constrains the set of  minimizers. As a result,
the  solution depends continuously on the input data. Instead of
minimizing (\ref{minimization}) over all functions, the minimization is
performed over a convex subset of functions, which results in a
stable reconstruction algorithm. The stabilization functional 
added to the minimization functional in (\ref{minimization}), 
with regularization parameter  $\alpha$, was
chosen \citep{Cherkaev:2001,Bonifasi:2009} as a quadratic functional
corresponding to Tikhonov regularization. We use this approach here
as well.






\subsection{Percolation theory}





Lattice and continuum percolation theories 
\citep{Stauffer-92,Torquato:RHM-02,Hornung:1997}
have been used to model a broad range of disordered
materials where the connectedness of one phase
dominates effective behavior. 
Consider the square ($d=2$) or cubic ($d=3$) network of bonds 
joining nearest neighbor sites on the integer lattice 
%${\mathbb{Z}}^d$.
${\mathbb{Z}}^d$.
%In the random pipe network,
We consider the problem of electrical transport 
through the network. 
The bonds are assigned electrical conductivities
$\sigma_0 > 0$ (open) or 0 (closed) with probabilities $p$ 
and $1-p$. 
Groups of connected open bonds 
are called open clusters, and the average cluster size
grows as $p$ increases.
The striking feature of this model is that 
there is a critical probability $p_c$, $0<p_c<1$,
called the {\it percolation threshold}, 
where an infinite cluster of open bonds first appears.
In $d=2$, $p_c = \frac{1}{2}$, and in $d=3$,
$p_c \approx 0.25$. 
Typical configurations for the $d=2$ square lattice
above and below the threshold
are shown in Figure \ref{percolation_theory}. 
Near the threshold, or transition point, the order parameters
of the system, such as the infinite cluster density,
the correlation length, and the effective conductivity
are characterized by universal critical exponents
which depend only on dimension, and not on the 
details of the lattice.





An important feature of phase transitions in statistical mechanics,
such as in the Ising model of a ferromagnet
\citep{Thompson-88},
is the collapse of a spectral gap as the critical point is approached
\citep{Baker-1990,Baker:PRL-990}.
Similarly, there is a gap in the support of the spectral measure
for $p$ below $p_c$
\citep{Clerc:AP-191,Day:JPCM-96,Golden:PRL-3935,%
Golden:JMP-5627,Jonckheere_Luck_JPA_1998},
which collapses as the threshold is approached. 










\section{Results}



\subsection{Numerical calculation of the spectral measure for bone}
%

In Section 2 
we discussed how the effective
properties of a composite such as bone may be expressed in terms of the
spectral measure $\mu$ of a self adjoint operator $\Gamma\chi$, where
$\chi$ is the indicator function for the bone and
$\Gamma=\nabla(-\Delta)^{-1}\nabla\cdot$ acts as a 
projector onto curl free fields. 
We observed that $\Gamma \chi$ is a 
self-adjoint operator with respect to the inner product 
with weight $\chi$. Here it is more convenient to use the 
standard $L^2$ inner product and consider the 
self-adjoint operator $\chi \Gamma \chi$ (the two approaches 
are easily seen to be equivalent).
In order to identify bone microstructures with spectral measures, 
we have set threshold levels and then discretized two 
SEM images of trabecular bone
from http://hansmalab.physics.ucsb.edu. 
The first image, shown in 
Figure \ref{bone_pics_spectral_measures}a, is from a
young (22-year old) healthy adult man. The second, shown in
Figure \ref{bone_pics_spectral_measures}b, is from 
an elderly, osteoporotic woman in her eighties.
Discretizations of these images are shown in 
Figure \ref{bone_pics_spectral_measures}c and 
\ref{bone_pics_spectral_measures}d.
We remark that the images in 
Figure \ref{bone_pics_spectral_measures}a and 
\ref{bone_pics_spectral_measures}b
are three dimensional, but our discretizations
are two dimensional representations (surface slices)
of the underlying images. Our subsequent calculations
of the spectral measure are then two dimensional
as well.


The discretizations are associated with networks 
of 1's (bone) and 0's (marrow), as in the 2-D square 
lattice percolation model, 
in order to calculate the corresponding spectral
measures. On this square lattice
the action of the differential operators $\nabla$ and $\nabla\cdot$ are 
defined in terms of the
forward and backward difference operators
\citep{Golden:CMP-467}. Furthermore, the action of $\chi$ is that of a
square diagonal matrix with ones and zeros along the diagonal
corresponding to type two bonds (bone) and type one bonds (marrow),
respectively. Therefore the action of the self adjoint operator
$\chi\Gamma\chi$ on this lattice is given by that of a real symmetric matrix,
and the spectral measure may be calculated directly from its
eigenvalues  $\{\lambda_i\}$ and eigenvectors $\{v_i\}$ via the formula
%\citep{Golden:CMP-473}  
%
\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}
%
where $\delta(\lambda)$ is the Dirac delta function and $e_0$ is a vector of
ones. The effective complex permittivity is given by 
%
\begin{equation}                    \label{Lattice_permittivity}
F(s)=\sum_i\frac{m_i}{s-\lambda_i}\;, \qquad
\epsilon^*=\epsilon_1(1-F(s)),
\end{equation}
%
%(see equation \eqref{eq:Fs}). 
where $s=1/(1-\epsilon_2/ \epsilon_1)$.


%The central object of the theory is the spectral function $\mu(\lambda)$
%which, for matrices, is defined pointwise at the eigenvalues with
%$m_i=\mu(\lambda_i)$. 




As the system size $N$ increases, the size of the matrix
$n$ increases proportionally by $N^2$, and the eigenvalues become
increasingly dense in the spectral interval [0,1]. 
For a random system averaged over many realizations 
(or using a large-enough spatial sample of a fixed realization), 
a high resolution histogram of the spectral
measure begins to resemble a smooth curve,
as in Figure \ref{RRN_realizations_spectral_measures}c and 
\ref{RRN_realizations_spectral_measures}d. 
The presence or absence of spectral
gaps at the endpoints of the spectral interval
and the details of how large the gap is or how large the spectral
values $m_i$ are at the collapse give important information pertaining to
the connectivity and electromagnetic transport properties of the system. 









We have approximated the 2-D geometric structure in pictures of sliced
healthy and osteoporotic trabecular bone by a $540\times540$ square
lattice bond network. The CPU time for solving the spectral problem for a
given square sample grows like $N^6$. We have optimized the
use of our algorithm by averaging over six coarser discretizations of
the original $540\times540$ discretizations. 
Each of the $90\times90$ coarser representations sample
every 6th data point of the original $540\times540$ system with a shifted
starting point, sampling every data point of the fine mesh $540\times540$
discretization. 
As in \citep{Jonckheere_Luck_JPA_1998}, where the spectral
measure was obtained by an indirect method, in our method
of calculating the spectral measure, we separate the calculation
of the delta function at the left endpoint $\lambda =0$
from the calculation of the spectral measure in the rest of the 
interval. In order to see the details of the spectral measure
throughout the interval, we avoid plotting the delta function
at $\lambda =0$.





The spectral measures for healthy and osteoporotic bone
are shown  
in Figure \ref{bone_pics_spectral_measures}e and 
\ref{bone_pics_spectral_measures}f, respectively.
The actual bone volume fractions in the discretized images are
0.5386 and 0.2431, respectively. For comparison,  
in Figure \ref{RRN_realizations_spectral_measures}a and 
\ref{RRN_realizations_spectral_measures}b 
we show realizations of the 
$30\times30$
percolation model with $p=0.5$
and $p=0.25$, respectively, along with direct calculations of the
spectral meaures (averaged over 5000 statistical samples)
in c and d.
The graphs show linearly connected peaks
of histograms with bin sizes on the order of $10^{-2}$ and $10^{-1}$ for
$p=0.5$ and $p=0.25$, respectively. The function constructed
in this way from the weights $m_i$ in the spectral measure
is called the {\it spectral function} $\mu(\lambda)$, and provides
a graphical representation of the measure.  
The apparent smoothness of
the spectral function graphs corresponding to the random bond networks
is due to the large number ($\sim10^6$) of eigenvalues and eigenvectors
calculated. In contrast, the jaggedness of the spectral function
graphs corresponding to the discretizations of trabecular bone is due
to the relatively small number ($\sim10^4$) of eigenvalues and
eigenvectors calculated. In order to reveal the detail of the spectral
functions with $p=0.5$ for the lattice model and for the young bone,
we have limited the plot domain to
$(\eta,1-\eta)$ with $\eta\sim10^{-15}$, eliminating some of the divergent
behavior in the graph. No such restriction was necessary for
the spectral functions with $p=0.25$ for either the lattice model
or the old bone. 

In Figure \ref{brine_connectivity} we show how the connectedness 
of the brine microstructure in sea ice changes dramatically 
with small changes in temperature. If trabecular connectedness
behaved similarly, an individual could develop osteoporosis
from a fever!







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%











\subsection{Numerical inversion for the spectral measure}










To simulate the effective complex permittivity of the sample, we
used data on the permittivity of bone and bone marrow in
\citep{Gabriel:1996a, Gabriel:1996b}, and calculated the values of
$F(s)$ for a set of points in the $s$-plane.  Some of the data 
in the s-plane used in the reconstruction are shown in Figure \ref{spoints}. 
The reconstructed spectral functions are shown in 
Figure \ref{spectral_function}. We use dashed
curves to indicate the results obtained for the old bone, and solid
curves to show the spectral functions corresponding to young bone. We
plot two curves in each case for reconstructions using the real or
imaginary parts of the integral operator $A$.

The porosity of the samples was estimated as the zeroth moment of the
spectral measure. We averaged the values over a range of the
regularization parameter $\alpha$ to avoid the problem of choosing
its appropriate value. We did not add any additional noise to the
data, however, the results obtained in \citep{Bonifasi:2009} indicate
that by using this approach, the values of the porosity can be stably
recovered from the effective properties, even in the presence of large noise.
In  \citep{Bonifasi:2009}, an L-curve method was used to find an optimal
value of the parameter $\alpha$, the regularized method used a Tikhonov stabilization functional.











\section{Discussion}





Properties of the spectral measure provide valuable information
about the structure of trabecular bone and the
degradation in the quality of connectedness as osteoporosis
progresses.  
Interestingly, we observe that the spectral function
for young, healthy bone with volume fraction $\phi \approx 0.54$
in Figure \ref{bone_pics_spectral_measures}e
closely resembles the spectral function for the lattice 
percolation model with $p=0.5$. Both geometries exhibit
long-range connectedness, and this is reflected in the
absence of a spectral gap in both cases. However, for the
osteoporotic bone with volume fraction $\phi \approx 0.24$, 
its spectral function shown in 
Figure \ref{bone_pics_spectral_measures}f
does not have an obvious gap as in the corresponding
lattice model with $p=0.25$. While the lattice model
is below its percolation threshold, the bone is still
connected over the scale of the image, even if the 
connections are tenuous or fragile,
and the gap is not present.
The tenuousness of
the connections is nevertheless manifested by a 
sudden decrease in the height of the spectral function 
from $m_i\sim10^{-6}$ for $1-\lambda_i\lesssim10^{-3}$ to
$m_i\sim10^{-35}$ for $1-\lambda_i\gtrsim10^{-3}$.  
Therefore, even though the gap has
effectively collapsed, the weakness of the connection is manifested in
the spectral function. Due to these tenuous connections the spectral
function of osteoporotic trabecular bone has characteristics of both
of the random networks shown in 
Figure \ref{RRN_realizations_spectral_measures}.




For our reconstructions of the spectral measure,
estimates of the bone volume fraction obtained by averaging 
the zeroth moment of the spectral
measures reconstructed with $\alpha = 0.1, 0.01, 0.001$  give
excellent agreement with the values calculated from the images. For
young and old bone samples, the true bone volume fraction 
was calculated as
0.5386  and  0.2431, and the reconstructed values were 0.5384 and
0.2418, respectively. Other features of the reconstructed 
spectral functions
for the young bone agree well with the spectral function in
Figure \ref{bone_pics_spectral_measures}e. The reconstructed
spectral function for the older bone, however, shows some 
discrepancy with 
Figure \ref{bone_pics_spectral_measures}f. While the zeroth
moments match up very closely, we believe that the different 
behavior near
$\lambda = 0$ may be due to the much larger bin size used
for the reconstructions. Such issues will be investigated further
in subsequent work.



We remark that the validity of our approach can be tested 
experimentally. The spectral measure for a bone
specimen can be calculated as we have here from
microstructural imagery. Then using this measure in an integral
representation for an effective material property,
such as the complex permittivity or electrical conductivity,
the predicted material property could be compared with 
experimental measurements of the same property. 
In this regard, while we have computed the spectral measure 
here only in two dimensions, the method extends easily
to three dimensions, although the numerical calculations
are significantly more intensive and time consuming.


In conclusion, we have introduced a new way of characterizing
the porous microstructure in bone. We have applied our
methods to specimens of trabecular bone, although the same
approach can be applied to the porous microstructure
of cortical bone as well. A particular strength of focusing on a 
spectral characterization of the microstructure as we have here
is that the spectral measure is precisely what enters into
calculations of material properties. On the other hand,
it is not clear at this stage how to obtain very specific,
small scale information about the microstructure, such as 
the sizes or shapes of the trabeculae from
the spectral measure. This represents a limitation of our
current approach, and an interesting area for future research.










\vspace{4ex}

\noindent 
{\bf Acknowledgments}\\ 
We are very grateful for the support provided by 
the Division of Mathematical Sciences and
the Arctic Natural Sciences Program
at the US National Science Foundation (NSF). 
Ben Murphy was supported by an
NSF VIGRE graduate fellowship, which
is gratefully acknowledged.
We would also like to thank David Dobson,
who wrote an early version of the 
program which computes the spectral measure,
and Maria-Grazia Ascenzi for her help in
providing numerous bone images as well as
with the manuscript. 


\vspace{4ex}




%\bibliographystyle{plain}



\bibliographystyle{plainnat}





\bibliography{golden,cmg,opp08,vert,bone,murphy,elena}



\eject













%%%%%%%%%  OLD Figures END



\noindent
Figure Captions


\vspace{4ex}





\noindent
Figure 1: 
The porous microstructure in the proximal femur of a
            52 year old female with normal bone density is shown
            on the left (photo courtesy of Maria-Grazia Ascenzi). 
        The brine inclusions in sea ice are shown on the right. 


\vspace{2ex}


\noindent
Figure 2: 
The two dimensional square bond lattice below
        its percolation threshold $p_c=1/2$ in a, 
        and above its threshold in b.
       Below $p_c$,
       there is no bulk transport, and above $p_c$
       the transport coefficient takes off with 
        power law behavior, as shown in c.

\vspace{2ex}




\noindent
Figure 3: 
SEM images of healthy and osteoporotic trabecular bone
            are shown in a and b, respectively, and their
            coarsened lattice discretizations are shown in c and d. 
            The graphs of the associated spectral functions
            are shown in e and f. The spectral
            functions have been normalized 
            by the uniform bin size so that the 
            area under each graph represents the 
            mass $\phi$ of the spectral measure, or the
            bone volume fraction.


\vspace{2ex}



\noindent
Figure 4: 
Realizations of the two dimensional lattice percolation 
            model are shown in a and b, and the corresponding
            spectral functions 
            (averaged over 5000 random realizations)
            are shown in c and d. In d, there is a spectral gap
            around $\lambda = 1$, indicating the lack of long-range
             order or connectedness. The gap collapses in c
            when the percolation threshold of $p=p_c=0.5$ has been reached,
            and the system exhibits long-range connectedness.
            Note the difference in vertical scale in the graphs 
            in c and d. 


\vspace{2ex}



\noindent
Figure 5: 
Thermal evolution of X-ray CT images of the brine phase 
within a lab-grown sea-ice single crystal 
with salinity of 9.3 ppt.  The (non-collocated)  
8 $\times$ 8 $\times$ 2 mm 
volumes illustrate a pronounced 
change in the micro-scale morphology and 
connectivity of the brine inclusions during warming.
As the sample warms beyond the percolation threshold
at a brine volume fraction $\phi$ of about 0.05,
the images show the development of long-range order, in the form of connected pathways
which facilitate transport. 




\vspace{2ex}



\noindent
Figure 6: 
Data points in the complex s-plane, where
         $s=1/(1-\epsilon_2/\epsilon_1)$, $\epsilon_1$
       is the complex permittivity of marrow, and 
        $\epsilon_2$ is the complex permittivity of bone.
       Both of these parameters depend on the frequency of the 
       applied electric field, and the values of $s$ form
       an arc in the complex plane as the frequency is varied
       over the range from 10 Hz to 100 GHz, although only the
       lower part of this range is shown here. 



\vspace{2ex}



\noindent
Figure 7: 
Reconstruction of the spectral function using 
            the Tikhonov regularization method. The dashed curves 
            denote reconstructions for the old bone, and the solid
            curves show the recovered spectral function 
            for the young bone. For these reconstructions, the 
           values have not been normalized by the bin size as 
           in Figures \ref{bone_pics_spectral_measures} 
          and \ref{RRN_realizations_spectral_measures}, so that
          the bone porosity is not obtained by the area under the
          curve, but by the sum of the reconstructed spectral data points. 
          Also, the bin size here is much larger than in 
          Figure \ref{bone_pics_spectral_measures}, 
         which explains the discrepancies in the details of the graphs.  






\eject





\begin{figure}          \centering
\includegraphics[height=4.5cm]{Figure_1.eps} 

%%%%%%%%%%%%  FIGURE 1

\caption{}

 \label{bone_sea_ice_microstructures}

\end{figure}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




\begin{figure}          \centering
\includegraphics[height=4.5cm] {Figure_2.eps}

%%%%%%%%%%%%  FIGURE 2

\caption{}

 \label{percolation_theory}

\end{figure}






\begin{figure}          \centering
\includegraphics {Figure_3.eps}

%%%%%%%%%%%%  FIGURE 3

\caption{}

 \label{bone_pics_spectral_measures}

\end{figure}



\begin{figure}          \centering
\includegraphics {Figure_4.eps}

%%%%%%%%%%%%  FIGURE 4

\caption{}

 \label{RRN_realizations_spectral_measures}

\end{figure}




\begin{figure}          \centering
\includegraphics {Figure_5.eps}

%%%%%%%%%%%%  FIGURE 5

\caption{}

 \label{brine_connectivity}

\end{figure}





\begin{figure}          \centering
\includegraphics [height=5cm] {Figure_6.eps}

%%%%%%%%%%%%  FIGURE 6

\caption{}

 \label{spoints}

\end{figure}








\begin{figure}          \centering
\includegraphics [height=7cm] 
   {Figure_7.eps}

%%%%%%%%%%%%  FIGURE 7

\caption{}

 \label{spectral_function}

\end{figure}






\end{document}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%







Let us briefly formulate the problem of finding the effective
or bulk conductivity of a two phase composite material
\citep{Golden:CMP-473,Golden:TGRS-1655,Gully:2007:Physica}.
For these considerations we ignore the air phase in sea ice.
Let the local conductivity 
$\sigma(x)$ be a spatially stationary random field
in $x \in \mathbb{R}^3$, for an appropriate probability space 
representing the set of realizations of the random medium.
While sea ice may exhibit significant variations in 
microstructural properties over its entire depth, 
relatively thin layers such as 5 or 10 centimeters
often display only small variations in these properties.
It is then reasonable to assume that in typical sample
sizes on this scale, the statistics describing the
brine microstructure in such a region are 
represented by a stationary random function throughout 
all of $\mathbb{R}^3$.
We assume $\sigma(x)$ 
takes the value $\sigma_1=\sigma_b$ 
in the brine phase, 
which depends on temperature $T$ via
\citep{Stogryn:IEEE-523}
\begin{align} \label{brineconductivity}
 \sigma_b=-T \exp(0.5193+0.08755T)\;\;\;
\Omega^{-1} \; \mbox{m}^{-1},   \\ \notag
T \geq -22.9^{\circ} \mbox{C}. \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\end{align}
In the ice phase, $\sigma(x)$ 
takes the value 
$\sigma_2=\sigma_i=1.1 \times 10^{-8} 
\; \Omega^{-1} \; \mbox{m}^{-1}$,
which is
the value at $T=-10^{\circ}$ C
\citep{Hobbs-74}. 
The value in ice is effectively 0 compared to the 
values of $\sigma_b$ in brine, which lie between
about 3 and 7 $\Omega^{-1} \; \mbox{m}^{-1}$.  
We write
%\begin{equation}                                 \label{eq2.1}
$\sigma = \sigma_1 \chi + \sigma_2 (1-\chi),$
%\end{equation}
where $\chi(x)$ is the characteristic function 
of the brine phase, with $\chi(x)=1$ for $x$ in the brine
and $\chi(x)=0$ for $x$ in the ice.  
The brine volume fraction $\phi$ is given by
$\phi = \langle \chi \rangle$, where
$\langle \cdot \rangle$ means
an ensemble average or spatial average over all of
$\mathbb{R}^3$, and depends on temperature $T$ ($^{\circ}$C)
and salinity $S$ parts per thousand (ppt) 
through the Frankenstein-Garner relation
\citep{Frankenstein:JG-943},
\begin{equation}       \label{Frankenstein_Garner}
 \phi=\frac{S}{1000}\left(\frac{49.185}{|T|}+0.532\right).
\end{equation}



Transport depends on the statistics of the brine inclusions, 
and in sea ice the inclusion cross-sectional areas $A$ 
(above a cut-off) have been observed to be lognormally distributed
\citep{Perovich:JGR-18327}. Then $Z=\ln{A}$ has a normal probability density
with mean $\alpha$ and variance $\delta^2$,
\begin{equation}                  \label{lognormal_distribution}
P(Z) = %(2 \pi \rho^2)^{-1/2}
     \frac{1}{\sqrt{2 \pi \delta^2}} \;
         e^{-(Z-\alpha)^2/2 \delta^2}.
\end{equation}
The mean brine inclusion horizontal 
cross-sectional radius increases with temperature $T$, or with
brine volume fraction $\phi$,
%\begin{equation}                  \label{mean_radius}
$r(\phi) = 7 \times 10^{-5} + (1.6 \times 10^{-4})\phi$ meters,
%\end{equation}
and also $\delta \approx 1$
\citep{Perovich:JGR-18327,Bock:2005,Golden:MM-801}.
The thermal evolution of brine phase geometry 
and connectivity is considerably more complicated than
is described by 
the lognormal distribution in (\ref{lognormal_distribution})
\citep{Golden:2007:GRL,Pringle:JGR:2009},
however, it does provide very useful information 
for understanding the transport properties of sea ice. 





Let $E(x)$ and $J(x)$ be the stationary
random electric and current fields satisfying
the constitutive law
%\begin{equation}                                 \label{eq2.2}
$    J (x)= \sigma(x) E (x) $
%\end{equation}
and the equations
\begin{equation}                                 \label{eq2.3}
\nabla \times E (x) = 0, \;\;\;\;\
\nabla \cdot J (x) = 0, 
\end{equation}
with
%\begin{equation}                                   \label{eq2.5}
$\langle E (x) \rangle = e_k,$
%\end{equation}
where $e_k$ is a unit vector in the $k^{th}$ direction for
some $k = 1, 2, 3$.
The effective conductivity tensor ${\bsig}^*$
is defined by \citep{Golden:CMP-473}
\begin{equation}                                    \label{eq2.6}
    \langle J \rangle= {\bsig}^* \langle E \rangle.
\end{equation}
For simplicity, we focus on one diagonal coefficient
$\sigma^* = \sigma^*_{kk}$. Due to the homogeneity of effective parameters,
$\sigma^*(\lambda \sigma_1, \lambda \sigma_2)=\lambda
\sigma^*(\sigma_1,\sigma_2)$, where $\lambda$ is any scalar,
$\sigma^*$ depends only on the ratio
$h = \sigma_1 / \sigma_2$, and we define
$m(h) =\sigma^* / \sigma_2\,$.
The two main properties of $m(h)$ are that it is analytic
off $(- \infty, 0]$ in the $h$--plane, and that it
maps the
upper half plane to the upper half plane,
so that it is an example of a Herglotz, or Stieltjes function.




An integral representation 
\citep{Bergman:PRC-377,Milton:APL-300,%
Golden:CMP-473,Golden:JMPS-333,Milton:CPAM-647}
for $m(h)$ which provides 
an important relationship between
microstructural information and the effective conductivity is
\begin{align}                                     	\label{eq2.8}
    \mathcal{F}(s)=1 - m(h)=\int_0^1 {\frac{d\mu(z)}{s-z}}\;, 
    \\ \notag  \qquad  s = 1/(1-h), \qquad  s\notin[0,1]\;,
\end{align}
where 
$\mathcal{F}(s)$
is analytic off $[0,1]$,
and $\mu$ is a positive measure on $[0,1]$.
Formula (\ref{eq2.8}) separates the
parameter $s$
from information about the mixture geometry 
contained in $\mu$, which is a spectral measure
of the operator $\Gamma\chi$, where
$\Gamma=\nabla(-\Delta)^{-1}\nabla \cdot$.
%is a projection operator onto curl-free fields.
Statistical assumptions about the random medium
(via the correlation functions) 
are incorporated through the moments $\mu_n$ of $\mu$.
For example, 
$\mu_0 = \int^1_0  d \mu (z) = 
\langle \chi \rangle = \phi,$ the porosity.













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




The key step in a rigorous method for obtaining bounds
on effective properties of composites, called the
{\it analytic continuation method} 
\cite{Bergman:PRC-377,Milton:APL-300,%
Golden:CMP-473,Golden:JMPS-333,Milton:CPAM-647},
is to obtain
an integral representation for $\sigma^*$.
This representation is quite powerful, and provides 
a useful and important relationship between
microstructural information and effective properties,
such as the bulk conductivity.
It is more convenient to work with 
the function 
%\begin{equation}                                    \label{eq2.7}
$\mathcal{F}(s) = 1 - m(h)$, where
$s = 1/(1-h)$,
%\end{equation}
which is analytic off $[0,1]$ in the $s$--plane.  
It was proven \cite{Golden:CMP-473} that $\mathcal{F}(s)$ has the 
representation
\begin{equation}                                     	\label{eq2.8}
    \mathcal{F}(s)=\int_0^1 {\frac{d\mu(z)}{s-z}}\;, \qquad  s\notin[0,1]\;,
\end{equation}
where $\mu$ is a positive measure on $[0,1]$.
Formula (\ref{eq2.8}) separates the
parameter information in $s$
from information about the mixture geometry 
contained in $\mu$, a spectral measure
of the operator $\Gamma\chi$, where
$\Gamma=\nabla(-\Delta)^{-1}\nabla \cdot$
is a projection operator onto curl-free fields.
Statistical assumptions about the 
geometry are incorporated into $\mu$ via its moments
\begin{equation}                                   \label{eq2.9}
\mu_n = \int^1_0 z^n d \mu (z) = 
(-1)^n \langle \chi 
[( \Gamma \chi )^n e_k ] \cdot e_k \rangle,
\end{equation}
which can be calculated from the
correlation functions of the random medium.
For example, the mass $\mu_0$ of $\mu$ is just the 
brine volume fraction,
\begin{equation}                                   \label{BVF}
\mu_0 = \int^1_0  d \mu (z) = 
\langle \chi \rangle = \phi.
\end{equation}

%For the complex elementary bounds it is assumed that 
%we know only $\mu_0 = p_1$, and if the medium
%is statistially isotropic we know
%$\mu_1 = p_1 p_2 / d$ as well.



Bounds on $\sigma^*$, or $\mathcal{F}(s)$, are obtained by fixing
$s$ in (\ref{eq2.8}), varying over admissible measures $\mu$ (or
admissible geometries), such as those that satisfy only $\mu_0=\phi$,
and finding the corresponding range of values of
$F(s)$ in the complex plane
\cite{Golden:CMP-473,Golden:JMPS-333}.
In general, $\sigma_1$
and $\sigma_2$ can be complex, although here we will consider
only real values for the component conductivities. 
The simplest bounds, where only the conductivities and 
volume fractions of the constituents are assumed known, 
are called the arithmetic (upper) and harmonic (lower) mean bounds,
\begin{equation}           \label{elementary_bounds}
\frac{1}{\frac{\phi}{\sigma_1}+\frac{1-\phi}{\sigma_2}} \;\;
\leq \;\; \sigma^* \;\; \leq \;\;
\phi \sigma_1 + (1-\phi) \sigma_2 \, .
\end{equation}
The bounds are optimal since laminates parallel to the field
realize the upper bound, 
while laminates perpendicular to the field realize the lower bound.
If the material is further assumed to be statistically
isotropic, i.e., $\sigma^*_{ik} = \sigma^* \delta_{ik}$,
then $\mu_1 = \phi (1-\phi) / d$, where $d$ is the dimensionality
where the material is isotropic (either $d=2$ or $d=3$).
Tighter inequalities known as the
Hashin-Shtrikman bounds \cite{Hashin:JAP-3125} are then obtained,
\begin{equation}					 \label{eq2.23}
\sigma_2  + \frac{\phi} { 
\frac{1}{\sigma_1 - \sigma_2} + \frac{1-\phi}{d
\sigma_2} } \;\; \leq \;\;  \sigma^* \;\; \leq \;\; 
\sigma_1 + \frac{1-\phi} { \frac{1}{\sigma_2
- \sigma_1} + \frac{\phi}{d \sigma_1} }\,,
\end{equation}
when $\sigma_1$ and $\sigma_2$ are real with 
$\sigma_1 \geq \sigma_2$.






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
















%Bounds on $\sigma^*$, or $\mathcal{F}(s)$, are obtained by fixing
%$s$ in (\ref{eq2.8}), varying over admissible measures $\mu$ (or
%admissible geometries), such as those that satisfy only $\mu_0=\phi$,
%and finding the corresponding range of values of
%$F(s)$ in the complex plane
%\cite{Golden:CMP-473,Golden:JMPS-333}.
%In general, $\sigma_1$
%and $\sigma_2$ can be complex, although here we will consider
%only real values for the component conductivities. 

Rigorous bounds on $\sigma^*$ can be obtained from
(\ref{eq2.8}) \citep{Golden:CMP-473,Golden:JMPS-333}.
The simplest bounds, where only the conductivities (assumed real) and 
volume fractions of the constituents are assumed known, 
are called the arithmetic (upper) and harmonic (lower) mean bounds,
\begin{equation}           \label{elementary_bounds}
\frac{1}{\frac{\phi}{\sigma_1}+\frac{1-\phi}{\sigma_2}} \;\;
\leq \;\; \sigma^* \;\; \leq \;\;
\phi \sigma_1 + (1-\phi) \sigma_2 \, .
\end{equation}
The bounds are optimal since laminates parallel to the field
realize the upper bound, 
while laminates perpendicular to the field realize the lower bound.
If the material is further assumed to be statistically
isotropic, i.e., $\sigma^*_{ik} = \sigma^* \delta_{ik}$,
then $\mu_1 = \phi (1-\phi) / d$, where $d$ is the dimensionality
where the material is isotropic (either $d=2$ or $d=3$).
Tighter inequalities known as the
Hashin-Shtrikman bounds \citep{Hashin:JAP-3125} are then obtained,
\begin{equation}					 \label{eq2.23}
\sigma_2  + \frac{\phi} { 
\frac{1}{\sigma_1 - \sigma_2} + \frac{1-\phi}{d
\sigma_2} } \;\; \leq \;\;  \sigma^* \;\; \leq \;\; 
\sigma_1 + \frac{1-\phi} { \frac{1}{\sigma_2
- \sigma_1} + \frac{\phi}{d \sigma_1} }\,,
\end{equation}
when $\sigma_1$ and $\sigma_2$ are real with 
$\sigma_1 \geq \sigma_2$.







\citep{Fritsch:JTB:2007}
‘Universal’ microstructural patterns in cortical and trabecular,
extracellular and extravascular bone materials: Micromechanics-based
prediction of anisotropic elasticity



\citep{Ural:JB:2006}
Cohesive finite element modeling of age-related toughness loss in
human cortical bone



\citep{Nalla:B:2004}
Effect of aging on the toughness of human cortical bone:
evaluation by R-curves


      

\citep{Hildebrand:JBMR:1999}
Direct Three-Dimensional Morphometric Analysis of
Human Cancellous Bone: Microstructural Data from
Spine, Femur, Iliac Crest, and Calcaneus





\citep{Oyen:JB:2008}
Composite bounds on the elastic modulus of bone


\citep{Coelho:JB:2009}
Numerical modeling of bone tissue adaptation — A hierarchical
approach for bone apparent density and trabecular structure


\citep{KnotheTate:BM:2009}
Bone as an inspiration for a novel class of mechanoactive materials

\citep{KnotheTate:JB:2008}
Idealization of pericellular fluid space geometry and dimension results in
a profound underprediction of nano-microscale stresses imparted by
fluid drag on osteocytes


\citep{Piekarski:N:1977} first fluid flow canaliculi



\citep{Issever:JBMR:2009} osteoporotic spine xray ct

\citep{Burghardt:JBMR:2010} Age- and Gender-Related 
Differences in the Geometric
Properties and Biomechanical Significance of
Intracortical Porosity



\citep{Goldenstein:ABE:2010}
Evaluation of the Presence of Bone Marrow in Cortical Porosity
in Postmenopausal Osteopenic Women



\citep{Hollister:JB:1994}
A Homogenization Sampling Procedure for
                  Calculating Trabecular Bone Effective
                  Stiffness and Tissue Level Stress


em props
\citep{Sierpowska:PMB:2003,Singh:CORR:1984,%
Saha:ITBE:1992,Williams:ABE:1996}

thermal conductivity
\citep{Davidson:MEP:2000,Biyikli:JBMR:1986}





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





Sea ice is also a porous composite whose properties depend
strongly on its microstructure. As sea water freezes,
brine is trapped in inclusions between pure ice grains,
forming a composite that also contains air and solid salts.
The polar sea ice packs play a fundamental
role in Earth's climate system, and are leading indicators
of climate change \citep{Thomas:2003:SI,Serreze:S-1533}. 
The brine microstructure in sea ice hosts extensive algal and 
bacterial communities which sustain life in the polar oceans
\citep{Thomas:2003:SI,Fritsen:S-782}. 
Fluid flow through the porous microstructure mediates 
a broad range of processes such as the growth and 
decay of seasonal ice, surface flooding and snow-ice formation, 
and the evolution of melt ponds and ice pack reflectance.
These processes are critical to our understanding of climate
change in the polar regions. The build-up of biomass in sea ice
depends on nutrient replenishment, which is
also constrained by fluid transport.
The fluid permeability of sea ice, in turn, depends
strongly on the volume fraction, or porosity, 
of the brine inclusions, as well as their connectedness.
These key microstructural characteristics depend 
strongly on temperature.

For example, below a porosity of about 5\%, sea ice is
effectively impermeable to fluid flow, yet is permeable for larger 
porosities \citep{Golden:S-2238,Golden:2007:GRL}. 
For a typical bulk salinity of 5 parts per thousand,
the critical porosity corresponds to a critical 
temperature of $-5^{\circ}$ C.
This critical behavior in sea ice is
due to a connectivity or percolation threshold 
in the brine microstructure occuring
at 5\% \citep{Golden:S-2238,Pringle:JGR:2009},
and is known as the {\it rule of fives}.


In recent work we have developed X-ray CT imaging methods
and pore structure analysis to characterize the thermal
evolution of the brine microstructure in sea ice
\citep{Golden:2007:GRL,Pringle:JGR:2009}.
We have brought to bear a suite of mathematical
techniques for predicting the temperature dependence
of the fluid permeability 
\citep{Golden:2007:GRL,Golden:MM-801,Zhu:AG-129,Golden:NAMS:2009}
and electrical conductivity of sea ice 
\citep{Zhu:PB:2010,Golden:2010:DSR}.
These techniques include 
percolation theory, hierarchical and network models, 
diffusion processes, and rigorous bounds.
We have also used analytic function theory to study
the effective complex permittivity of sea ice 
for remote sensing applications
\citep{Golden:JGR-13699,Golden:IMA-97,Golden:TGRS-1655,%
Golden:TGRS-1675,Gully:2007:Physica}, and inverted
permittivity data to obtain bounds on the brine porosity
of sea ice \citep{Cherkaeva:WRM-437}.


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


Forward bounds on $\e^*$ in the complex plane
are obtained by fixing
$s$ in (\ref{eq2.8}), varying over admissible measures $\mu$ 
(or admissible geometries), 
such as those that satisfy only $\mu_0=\phi$,
and finding the corresponding range of values of
$F(s)$ in the complex plane. As more information
is known about the microgeometry, the bounds become 
tighter, and a series of nested, lens-shaped 
regions is obtained. In the case of real
component parameters, the classical arithmetic and 
harmonic mean bounds are obtained if one only assumes
knowledge of $\phi$. The Hashin-Shtrikman bounds
are obtained if one further assumes statistical isotropy
\citep{Hashin:JAP-3125}. The complex bounds, as well as the
real case, are 
covered thoroughly in 
\citep{Golden:IMA-97,Golden:TGRS-1655,Gully:2007:Physica}.



The objective of inverse bounds, which is one approach
within {\it inverse homogenization},
is to use data about 
the electromagnetic response of a composite such
sea ice or bone to recover 
information about its structural parameters, in particular 
the brine volume fraction $\phi$, or bone volume fraction.  
The inverse method developed in 
\citep{Cherkaeva:WRM-437,Golden:TGRS-1675,%
McPhedran:1982,McPhedran:1990} 
yields intervals of uncertainty for the brine or bone
volume fraction $\phi$.  
Given an observed value of the complex 
permittivity $\e^*$, $\phi$ 
is increased until the 
value of $\e^*$ touches one boundary of 
a bounding region,
and is then decreased until the value touches the 
other boundary.
This procedure gives a range of 
values $\phi^{\ell}_1 \leq \phi \leq \phi^u_1$,
which can be made smaller by inverting tighter
forward bounds.
Below we will use a more sophisticated inversion
scheme over frequency
to reconstruct the spectral measure, and obtain
accurate estimates of bone volume
fraction $\phi$.


Examples include rocks 
\citep{Bourbie:JGR-11524,Broadbent:PCPS-629}, 
semiconductors \citep{Efros-84}, thin films 
\citep{Davis:OC-70}, glacial ice \citep{Enting:1985:LSM}, 
polycrystalline metals \citep{Chen:PRL:2007}, 
radar absorbing coatings \citep{Kusy:N-58}, and 
carbon nanotube composites \citep{Kyrylyuk:PNAS:2008}.  






The infinite cluster density 
$P_{\infty}(p)$ is defined as the probability that the origin 
(or any point, by translation invariance)
is contained in the infinte cluster.
The two-point function $\tau(x,y)$ measures the probability
that two points $x$ and $y$ belong to the same open
cluster. As $|x-y| \rightarrow \infty$, 
$\tau(x,y) \sim exp(-|x-y|/\xi)$, where $\xi$ is called the 
correlation length, which is a measure of the scale of connectedness. 
The correlation length diverges as $p$ increases toward $p_c$
and the infinite cluster forms.
In the neighborhood of $p_c$, $P_{\infty}(p)$ and $\xi(p)$
exhibit power laws, with 
%\begin{equation}
$P_{\infty}(p) \sim (p-p_c)^{\beta}$, $ p \rightarrow p_c^+$,
and
$\xi(p) \sim (p_c - p)^{-\nu}$, $p \rightarrow p_c^-$,
%\end{equation}
where $\beta$ and $\nu$ are {\it universal} critical exponents which
are believed to depend only on dimension, and not on the details
of the lattice. In $d=3$, $\beta \approx 0.41$ and $\nu \approx 0.88$
\citep{Stauffer-92}.








Let $\sigma^*(p)$ be the effective conductivity 
of this random resistor network
in the vertical direction, as defined by Kirchoff's laws
in the infinite volume limit
\citep{Golden:PRL-2923}.  
For $p < p_c$, $\sigma^*(p)=0$, as shown in 
Figure \ref{percolation_theory} (c). For $p>p_c$ near the threshold,  
$\sigma^*(p)$ exhibits power law behavior,
\begin{equation}                    \label{critical_conductivity}
\sigma^*(p) \sim \sigma_0(p-p_c)^t\;,\;\;\;\; p \rightarrow p_c^+,
\end{equation}
where $t$ is the conductivity critical exponent.
%There was debate concerning the value of $t$ in three dimensions  
%over roughly three decades, 
%with a wide range of published 
For lattices, %$p_c$ depends on the type of lattice, while 
$t$ is believed to be universal, 
depending only on $d$.
In $d=2$, $t \approx 1.3$,
and in $d=3$,
$t \approx 2.0$ \citep{Stauffer-92,Torquato:RHM-02}.  
There is also a rigorous 
bound \citep{Golden:PRL-2923,Golden:CMP-467} that $1 \leq t \leq 2$
in $d=2$ and $d=3$.







