Stieltjes functions and spectral analysis in the physics of sea ice

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, 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 how to use information on smaller-scale structures 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. Motivated by sea ice remote sensing, the analytic continuation method for obtaining rigorous bounds on the homogenized coefficients of two-phase composites was applied to the complex permittivity of sea ice, which is a Stieltjes function of the ratio of the 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 homogenization, 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 and providing new insights into the physics of sea ice.

1 Introduction

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 (Stroeve et al., 2007, 2012; Maslanik et al., 2007; Notz and Community, 2020; Notz and Stroeve, 2016). 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 (Turner et al., 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 (Golden et al., 2020).

Figure 1 Sea ice as a multiscale composite material. From left to right: millimeter-scale brine inclusions that form the porous microstructure of sea ice (Golden et al., 2007); centimeter-scale polycrystalline structure of sea ice (Arcone et al., 1986); melt ponds on Arctic sea ice in late spring and summer (Don Perovich) turning the surface into a two-phase composite of ice and melt water; the sea ice pack as a granular composite viewed from space (NASA), with “grains” ranging in horizontal extent from meters to tens of kilometers; and the Arctic Ocean viewed from space (NASA).

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 10 orders of magnitude, from the sub-millimeter scale to thousands of kilometers, as indicated in Fig. 1. Modeling the macroscopic behavior of sea ice on scales appropriate for climate models and process studies 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 resolution to tens of kilometers and typically require clever approximations and parameterizations to incorporate the basic physics of sea ice (Golden et al., 2020; Golden, 2015, 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 (Thomas and Dieckmann, 2003; Reimer et al., 2022).

We observe that this central problem of finding the effective properties of sea ice 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 (Thompson, 1988; Christensen and Moloney, 2005). Moreover, mathematical homogenization theory for differential equations with random coefficients similarly seeks to find the large-scale effective behavior from some knowledge about the local coefficients (Milton, 2002; Torquato, 2002; Bensoussan et al., 1978; Papanicolaou and Varadhan, 1982; Kozlov, 1989). 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) (Bergman, 1980; Milton, 1980; Golden and Papanicolaou, 1983; Golden, 1997b; Milton, 2002), in particular, yields powerful integral representations for the effective or homogenized transport coefficients of two-component (Golden and Papanicolaou, 1983) or multicomponent (Golden and Papanicolaou, 1985; Golden, 1986) 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 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 spectral measure of G. 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 (Murphy et al., 2015).

The Stieltjes or Herglotz structure of the effective parameters and their integral representations can be used to find rigorous bounds on the homogenized transport coefficients (Bergman, 1980; Milton, 1980; Golden and Papanicolaou, 1983; Golden, 1986; Baker and Graves-Morris, 1996; Milton, 2002; Gully et al., 2015; Cherkaev, 2019). 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 (Golden, 1995; Golden et al., 1998b, c). For example, the mass of the spectral measure is the brine volume fraction. If this is known, then one can obtain 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 a matrix-particle structure, e.g., separated brine inclusions in a pure ice host (Bruno, 1991; Golden, 1997b). 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 (Cherkaev and Golden, 1998; Cherkaev, 2001), one uses measurements of bulk electromagnetic behavior, such as the effective complex permittivity, to obtain information about the spectral measure (Cherkaev, 2001) and bounds on the microstructural characteristics, such as the brine volume fraction, connectivity (Cherkaev and Tripp, 1996; Cherkaev and Golden, 1998; Golden et al., 1998b; Gully et al., 2007; Orum et al., 2012; Cherkaev and Bonifasi-Lista, 2011), and crystal orientation (Gully et al., 2015). The microscale structure, which determines the spectral measure and the homogenized coefficient, is thus linked to the macroscopic behavior via the operator G and its spectral characteristics, and 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 (Golden and Papanicolaou, 1985; Golden, 1986).

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 (McLaughlin et al., 1985; Biferale et al., 1995; Fannjiang and Papanicolaou, 1994, 1997; Pavliotis, 2002; Majda and Kramer, 1999; Majda and Souganidis, 1994; Xin, 2009). 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 (1989, 1991) found a Stieltjes integral representation for the effective diffusivity as a function of the Péclet number for diffusion in an incompressible velocity field. Based on the approach in Golden and Papanicolaou (1983), 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 Murphy et al. (2017b, 2020a) we proved novel versions of the Stieltjes 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 (Bergman, 1982; Golden, 1986) for both the time-dependent and time-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 Gully et al. (2015), based on a resolvent formula for the electric field and earlier observations in Milton (1981), Bergman and Stroud (1992), and Milton (2002). 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 (Gully et al., 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 (Weeks and Ackley, 1982; Petrich and Eicken, 2009; Untersteiner, 1986). For example, while sea ice grown in quiescent conditions tends to have rather large-grained columnar structure, when grown in more turbulent or wavy conditions, it typically has a fine-grained granular structure. These distinctly different ice types have quite different fluid flow properties (Golden et al., 1998a). Additionally, when there is a well-defined current direction during formation, crystal orientations tend to be statistically anisotropic within the horizontal plane (Weeks and Gow, 1980). This can significantly affect the sea ice radar signature used in measurements of sea ice thickness and other properties used to validate climate models (Golden and Ackley, 1981).

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 (Waseda et al., 2018). 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 (Strong and Rigor, 2013). These recent changes can have complex implications for both sea ice formation and melting (Li et al., 2021). 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 model the surface layer of the ice-covered ocean as a continuum with effective properties (Bates and Shapiro, 1980; Keller, 1998; Wang and Shen, 2010; Mosig et al., 2015). Recently, 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 (Sampson, 2017). 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 sea ice concentration field. If the measure's mass is known, rigorous bounds can be obtained on the complex viscoelastic parameter. Previously, this effective parameter had only been fitted to wave data.

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 (Murphy et al., 2015) and studying the properties of computed spectral measures for a broad range of sea ice and other microstructures, like human bone (Golden et al., 2011; Bonifasi-Lista and Cherkaev, 2008, 2009; Bonifasi-Lista et al., 2009; Cherkaev and Bonifasi-Lista, 2011), 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) (Anderson, 1958; Evers and Mirlin, 2008), which provides a powerful theoretical framework for understanding 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 (Shklovskii et al., 1993; Kravtsov and Muttalib, 1997). For small disorder, the wave functions are extended and overlap, giving rise to correlated Wigner–Dyson (WD) statistics (Shklovskii et al., 1993; Kravtsov and Muttalib, 1997) with strong level repulsion (Guhr et al., 1998). In work on 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), we found an Anderson transition in spectral characteristics as the microstructure developed long-range order in the approach to a percolation threshold (Murphy et al., 2017a). 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, 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 (Guhr et al., 1998). In Morison et al. (2022) a novel class of two-phase media was introduced – twisted bilayer composites – based on Moiré 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 (Murphy and Golden, 2012; Murphy et al., 2015; Gully et al., 2015; Murphy et al., 2017a, b, 2020a; Golden et al., 2020; Golden, 2015, 2009, 1997b; Golden et al., 1998b, c). We also mention a significant recent advance in obtaining a Stieltjes integral representation for the fluid permeability of a porous medium (Bi et al., 2023) and an excellent, recent review of Stieltjes integrals in materials science (Luger and Ou, 2022; Ou and Luger, 2022). The permeability result is relevant for sea ice modeling (Golden et al., 1998a, 2007; Golden, 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 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. 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.

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. The interplay between homogenization techniques like the analytic continuation method here and models of phase transitions in statistical physics (Banwell et al., 2023) is particularly interesting across the full range of scales. From the millimeter-scale brine inclusions (Golden et al., 1998a, 2007) to meter-scale melt ponds (Ma et al., 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. 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, the theory holds in the quasistatic regime where the wavelength in the medium is assumed to be much longer than the microstructural scale. Typically, electromagnetic 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 (Morison et al., 2022), biomedical engineering of human bone (Golden et al., 2011; Bonifasi-Lista and Cherkaev, 2008, 2009; Bonifasi-Lista et al., 2009; Cherkaev and Bonifasi-Lista, 2011), and our understanding of the physical properties of general polycrystalline media in geosciences and other fields (Gully et al., 2015; Cherkaev, 2019).

2 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 (Broadbent and Hammersley, 1957; Stauffer and Aharony, 1992; Grimmett, 1989; Bunde and Havlin, 1991) focuses on connectedness in disordered and inhomogeneous media and has provided the theoretical framework for describing the behavior of fluid flow through sea ice (Golden et al., 1998a, 2007; Golden, 2009).

Consider the d-dimensional integer lattice ℤ^d, and the square or cubic network of bonds joining nearest-neighbor lattice sites. In the percolation model (Broadbent and Hammersley, 1957; Stauffer and Aharony, 1992; Grimmett, 1989; Bunde and Havlin, 1991), we assign to each bond a conductivity σ₀>0 with probability p, meaning it is open (black), and with probability 1−p we assign σ₀=0, meaning it is closed. Two examples of lattice configurations are shown in Fig. 2 with $p=\mathrm{1}/\mathrm{3}$ in (a) and $p=\mathrm{2}/\mathrm{3}$ in (b). Groups of connected open bonds are called open clusters. In this model there is a critical probability p_c, $\mathrm{0}<p_{\mathrm{c}}<\mathrm{1}$, the percolation threshold, at which the average cluster size diverges and an infinite cluster appears. For the d=2 bond lattice, $p_{\mathrm{c}}=\mathrm{1}/\mathrm{2}$. For p<p_c the infinite cluster density P_∞(p)=0, while for p>p_c, P_∞(p)>0 and near the threshold, $P_{\mathrm{\infty }}\left(p\right)\sim (p-p_{\mathrm{c}}{)}^{\mathit{\beta }}$ as $p\to p_{\mathrm{c}}^{+}$, where β is a universal critical exponent. It depends only on dimension and not on the details of the lattice. Let $\mathit{x},\mathit{y}\in {\mathbb{Z}}^d$ and τ(x,y) be the probability that x and y belong to the same open cluster. Then for p<p_c, $\mathit{\tau }(\mathit{x},\mathit{y})\sim e^{-|\mathit{x}-\mathit{y}|/\mathit{\xi }\left(p\right)}$, and the correlation length $\mathit{\xi }\left(p\right)\sim (p_{\mathrm{c}}-p{)}^{-\mathit{\nu }}$ diverges with a universal critical exponent ν as $p\to p_{\mathrm{c}}^{-}$, as shown in Fig. 2c.

The effective conductivity σ^*(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_∞(p) since there are no infinite pathways, as shown in Fig. 2e. For p>p_c, ${\mathit{\sigma }}^{*}\left(p\right)>\mathrm{0}$, and near p_c, ${\mathit{\sigma }}^{*}\left(p\right)\sim {\mathit{\sigma }}_{\mathrm{0}}(p-p_{\mathrm{c}}{)}^t,p\to p_{\mathrm{c}}^{+},$ where t is the conductivity critical exponent, with $\mathrm{1}\le t\le \mathrm{2}$ in $d=\mathrm{2},\mathrm{3}$ (Golden, 1990, 1992, 1997a), and numerical values t≈1.3 in d=2 and t≈2.0 in d=3 (Stauffer and Aharony, 1992). Consider a random pipe network with fluid permeability k^*(p) exhibiting similar behavior $k^{*}\left(p\right)\sim k_{\mathrm{0}}(p-p_{\mathrm{c}}{)}^e,$ where e is the permeability critical exponent, with e=t (Chayes and Chayes, 1986; Sahimi, 1995; Golden, 1997a). 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≠t (Halperin et al., 1985; Feng et al., 1987; Stauffer and Aharony, 1992; Sahimi, 1994; Kerstein, 1983).

Figure 2 The two-dimensional square lattice percolation model below its percolation threshold of $p_{\mathrm{c}}=\mathrm{1}/\mathrm{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).

3 Analytic continuation for two-phase composites

We now describe the analytic continuation method (ACM) for studying the effective properties of composites (Bergman, 1980; Milton, 1980; Golden and Papanicolaou, 1983; Golden, 1997b). 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 (Cherkaev, 2001; Cherkaev and Bonifasi-Lista, 2011; Cherkaev, 2019; Ou and Cherkaev, 2006; Ou, 2012; Kantor and Bergman, 1982; Avellaneda and Majda, 1989, 1991; Murphy et al., 2020b, 2017b). In the works of Golden (Golden, 1995, 1997b, 2015, 2009; Golden et al., 1998b, c, 2020), rigorous bounds on the complex permittivity of sea ice were found.

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 local permittivity tensor ϵ(x,ω), a spatially stationary random field in x∈ℝ^d and ω∈Ω, where Ω is the set of realizations of the medium. We consider a two-phase locally isotropic medium, where the components ϵ_jk, $j,k=\mathrm{1},..,d$, of ϵ satisfy

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{ϵ}}_{jk}(\mathit{x},\mathit{\omega })=\mathit{ϵ}(\mathit{x},\mathit{\omega }){\mathit{\delta }}_{jk},\end{array}$$

where d is dimension, δ_jk is the Kronecker delta, and

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{ϵ}(\mathit{x},\mathit{\omega })={\mathit{ϵ}}_{\mathrm{1}}{\mathit{\chi }}_{\mathrm{1}}(\mathit{x},\mathit{\omega })+{\mathit{ϵ}}_{\mathrm{2}}{\mathit{\chi }}_{\mathrm{2}}(\mathit{x},\mathit{\omega }).\end{array}$$

Later, we will consider a polycrystalline medium where ϵ is a non-trivial symmetric matrix. Here, χ_i(x,ω) is the characteristic function of medium $i=\mathrm{1},\mathrm{2}$, equaling 1 for ω∈Ω with medium i at x and 0 otherwise, with ${\mathit{\chi }}_{\mathrm{1}}+{\mathit{\chi }}_{\mathrm{2}}=\mathrm{1}$. The random electric and displacement fields E(x,ω) and D(x,ω) satisfy

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathbf{\nabla }\times \mathit{E}=\mathrm{0},\mathbf{\nabla }\cdot \mathit{D}=\mathrm{0},\mathit{D}=\mathit{ϵ}\mathit{E}.\end{array}$$

A variational problem (Golden and Papanicolaou, 1983) establishes that E can be written as $\mathit{E}={\mathit{E}}_f+{\mathit{E}}_{\mathrm{0}}$ satisfying

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathit{E}={\mathit{E}}_f+{\mathit{E}}_{\mathrm{0}},\mathbf{\nabla }\times {\mathit{E}}_f=\mathrm{0},\langle \mathit{D}\cdot {\mathit{E}}_f\rangle =\mathrm{0},\langle \mathit{E}\rangle ={\mathit{E}}_{\mathrm{0}}.\end{array}$$

In simpler terms, this says that curl-free and divergence-free fields are orthogonal (Helmholtz's theorem).

The effective permittivity tensor ϵ^* is defined as $\langle \mathit{D}\rangle ={\mathit{ϵ}}^{*}\langle \mathit{E}\rangle$, where 〈⋅〉 is ensemble averaging over Ω or, by an ergodic theorem, spatial average over all of ℝ^d (Golden and Papanicolaou, 1983). We prescribe that E₀ has direction e_k, the kth direction unit vector, and focus on the diagonal coefficient ${\mathit{ϵ}}^{*}={\mathit{ϵ}}_{kk}^{*}$, with ${\mathit{ϵ}}^{*}=\langle \mathit{ϵ}\mathit{E}\cdot {\mathit{e}}_k\rangle$. The key step of the method is to obtain the following Stieltjes integral representation for ϵ^* (Bergman, 1978; Milton, 1980; Golden and Papanicolaou, 1983; Milton, 2002):

$$\begin{array}{}\text{(5)}& F\left(s\right)=\mathrm{1}-{\displaystyle \frac{{\mathit{ϵ}}^{*}}{{\mathit{ϵ}}_{\mathrm{2}}}}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{\mathrm{d}\mathit{\mu }\left(\mathit{\lambda }\right)}{s-\mathit{\lambda }}},s={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{ϵ}}_{\mathrm{1}}/{\mathit{ϵ}}_{\mathrm{2}}}},\end{array}$$

where μ is a positive Stieltjes measure with support in [0,1], and F plays the role of a (negative) electric susceptibility (Bergman, 1978). In the variable $s=\mathrm{1}/(\mathrm{1}-h)$, with $h={\mathit{ϵ}}_{\mathrm{1}}/{\mathit{ϵ}}_{\mathrm{2}}$, F(s) is a Stieltjes function (Golden, 1997c; Cherkaev, 2001; Murphy and Golden, 2012). This representation arises from a resolvent formula for the electric field (in medium 1) (Murphy et al., 2015; Cherkaev, 2001; Golden and Papanicolaou, 1983):

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\chi }}_{\mathrm{1}}\mathit{E}=s(sI-G{)}^{-\mathrm{1}}{\mathit{\chi }}_{\mathrm{1}}{\mathit{e}}_k,G={\mathit{\chi }}_{\mathrm{1}}\mathrm{\Gamma }{\mathit{\chi }}_{\mathrm{1}},\end{array}$$

yielding $F\left(s\right)=\langle \left[\right(sI-G{)}^{-\mathrm{1}}{\mathit{\chi }}_{\mathrm{1}}{\mathit{e}}_k]\cdot {\mathit{e}}_k\rangle$, where $\mathrm{\Gamma }=-\mathbf{\nabla }(-\mathrm{\Delta }{)}^{-\mathrm{1}}\mathbf{\nabla }\cdot$ is a projection onto the range of the gradient operator ∇ and I is the identity operator. Equation (5) is the spectral representation of the resolvent formula in Eq. (6), and μ is a spectral measure of the self-adjoint operator G=χ₁Γχ₁ on L²(Ω,P).

A critical feature of Eq. (5) is that the component parameters in s are separated from the geometrical information in μ. Information about the geometry enters through the moments

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{\mu }}_n=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\mathit{\lambda }}^n\mathrm{d}\mathit{\mu }\left(\mathit{\lambda }\right)=\langle G^n{\mathit{\chi }}_{\mathrm{1}}{\mathit{e}}_k\cdot {\mathit{\chi }}_{\mathrm{1}}{\mathit{e}}_k\rangle .\end{array}$$

Then μ₀=ϕ, where ϕ is the volume or area fraction of phase 1, such as the brine volume fraction, the open-water area fraction, or melt pond coverage, and ${\mathit{\mu }}_{\mathrm{1}}=\mathit{\varphi }(\mathrm{1}-\mathit{\varphi })/d$ if the material is statistically isotropic. In general, μ_n depends on the (n+1)-point correlation function of the medium. This integral representation yields rigorous forward bounds for the effective parameters of composites, given partial information on the microgeometry via the μ_n (Bergman, 1980; Milton, 1980; Golden and Papanicolaou, 1983; Bergman, 1982). The integral representations can also yield 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 (McPhedran et al., 1982; McPhedran and Milton, 1990; Cherkaev and Tripp, 1996; Cherkaev and Golden, 1998; Cherkaev, 2001; Zhang and Cherkaev, 2009; Bonifasi-Lista and Cherkaev, 2009; Cherkaev and Bonifasi-Lista, 2011; Day and Thorpe, 1999; Golden et al., 2011; Cherkaev, 2020). See Sect. 5 for more details.

3.1 Spectral measure computations for two-phase composites

Computing the spectral measure μ for a given 2D composite microstructure first involves discretizing a two-phase image of the composite into a square lattice filled with 1s and 0s corresponding to the two phases. On this square lattice, the action of the differential operators ∇ and ∇⋅ is defined in terms of forward- and backward-difference operators (Golden, 1992; Murphy et al., 2015). Then the key operator χ₁Γχ₁, which depends on the geometry of the network via χ₁, becomes a real-symmetric matrix M (Murphy et al., 2015). Here Γ is a (non-random) projection matrix which depends only on the lattice topology and boundary conditions, and χ₁ is a diagonal (random) projection matrix which determines the geometry and component connectivity of the composite medium (Gully et al., 2015; Murphy et al., 2015). Another spectral approach to finding effective properties based on analytic continuation relies on computation of the electromagnetic eigenstates of individual inclusions (Bergman et al., 2020; Bergman, 2022).

The powerful integral representation in Eq. (5) 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 ∇, ∇×, and ∇⋅ (Huang et al., 2019) and will be published elsewhere. The integral representation in Eq. (8) is established in Theorem 2.1 of Murphy et al. (2015).

Theorem 1.

For each ω∈Ω, let M(ω)=W(ω)Λ(ω) W^T(ω) be the eigenvalue decomposition of the real-symmetric matrix M(ω)=χ₁(ω) Γ χ₁(ω). Here, the columns of the matrix W(ω) consist of the orthonormal eigenvectors w_i(ω), $i\mathrm{=}\mathrm{1}\mathrm{,}\mathrm{\dots }\mathrm{,}N$, of M(ω), and the diagonal matrix $\mathbf{\Lambda }\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\text{diag}\mathrm{\left(}{\mathit{\lambda }}_{\mathrm{1}}\mathrm{\right(}\mathit{\omega }\mathrm{)}\mathrm{,}\mathrm{\dots }\mathrm{,}{\mathit{\lambda }}_N\mathrm{(}\mathit{\omega }\mathrm{\left)}\mathrm{\right)}$ involves its eigenvalues λ_i(ω). Denote ${\mathbf{Q}}_i\mathrm{=}{\mathit{w}}_i{\mathit{w}}_i^T$ the projection matrix onto the eigenspace spanned by w_i and denote ${\mathit{\delta }}_{{\mathit{\lambda }}_i}\mathrm{\left(}\mathrm{d}\mathit{\lambda }\mathrm{\right)}$ the Dirac δ measure centered at λ_i. The electric field E(ω) satisfies $\mathit{E}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}{\mathit{E}}_{\mathrm{0}}\mathrm{+}{\mathit{E}}_f\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$, with E₀=〈E(ω)〉 and ΓE(ω)=E_f(ω), and χ₁E satisfies the resolvent formula in Eq. (6) with G replaced by M. The effective complex permittivity tensor ϵ^* has components ${\mathit{ϵ}}_{jk}^{\mathrm{*}}$, $j\mathrm{,}k\mathrm{=}\mathrm{1}\mathrm{,}\mathrm{\dots }\mathrm{,}d$, which satisfy

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\mathit{ϵ}}_{jk}^{*}={\mathit{ϵ}}_{\mathrm{2}}({\mathit{\delta }}_{jk}-F_{jk}(s\left)\right),F_{jk}\left(s\right)=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{\mathrm{d}{\mathit{\mu }}_{jk}\left(\mathit{\lambda }\right)}{s-\mathit{\lambda }}},\\ & \mathrm{d}{\mathit{\mu }}_{jk}\left(\mathit{\lambda }\right)=\sum _{i=\mathrm{1}}^N\langle {\mathit{\delta }}_{{\mathit{\lambda }}_i}\left(\mathrm{d}\mathit{\lambda }\right){\mathit{\chi }}_{\mathrm{1}}{\mathbf{Q}}_i{\widehat{\mathit{e}}}_j\cdot {\widehat{\mathit{e}}}_k\rangle .\end{array}\end{array}$$

Figure 3 Electric fields and spectral measures for sea ice brine microstructure. (a) Electric field values in log-10 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 μ(λ) (histogram representations) of the spectral measure μ in (c) and (d), which display spectral weights $m_i^{\mathrm{2}}$ versus the associated eigenvalues λ_i of the matrix M and index i, respectively, where ϕ is the brine area fraction of the images in (a). The vertical bars in (d) delineate the δ functions at the spectral endpoints $\mathit{\lambda }=\mathrm{0},\mathrm{1}$ seen in (c) from the rest of the spectrum, where ${\mathit{\lambda }}_i\lesssim {\mathrm{10}}^{-\mathrm{14}}$ and $\mathrm{1}-{\mathit{\lambda }}_i\lesssim {\mathrm{10}}^{-\mathrm{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 δ function at λ=0 switching from numerically zero, with the $m_i^{\mathrm{2}}\lesssim {\mathrm{10}}^{-\mathrm{30}}$, to $m_i^{\mathrm{2}}\sim \mathrm{1}$, giving rise to the “hot spots” in E via Eq. (9). The electrical permittivity is taken to be ${\mathit{\epsilon }}_{\mathrm{1}}=\mathrm{63.3}+i\mathrm{1930}$ for brine and ε₂=3.06 for ice (Backstrom and Eicken, 2006). E₀ is taken to be vertically oriented.

Figure 4 Temperature gradient fields and spectral functions for sea ice melt pond microstructure. (a) Temperature gradient field values in log-10 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 μ(λ) displayed below in (b). Here the gradient of temperature T, −∇T, plays the role of E_f in Eqs. (3) and (4) (Milton, 2002) (${\mathit{E}}_f=-\mathbf{\nabla }\mathit{\phi }$ for some electrical potential φ); the component thermal conductivities κ_i, $i=\mathrm{1},\mathrm{2}$, play the role of the complex electrical permittivities ϵ_i; and the heat current Q plays the role of the displacement field D; here, ϕ is the melt pond area fraction for the images in (a). Analogous to Fig. 3, as the fluid phase becomes increasingly connected on macroscopic length scales, a buildup of spectral measure mass at λ=0 shown in (b) leads to the formation of a δ function at λ=0, with corresponding switching in the values of the $m_i^{\mathrm{2}}$ from numerically zero, with the $m_i^{\mathrm{2}}\lesssim {\mathrm{10}}^{-\mathrm{30}}$ for the left an middle figure panels, to $m_i^{\mathrm{2}}\sim \mathrm{1}$ for the rightmost panel. The δ function at λ=1 is also analogous to that in Fig. 3. Like E₀ in Fig. 3, we take the average thermal gradient to be vertically oriented. The thermal conductivity is taken to be κ₁=0.5606 W m⁻¹ K⁻¹ for melt ponds and κ₂=0.3073 W m⁻¹ K⁻¹ for the surrounding snow (Yen, 1981).

Figure 5 Temperature gradient fields and spectral functions for Arctic pack ice microstructure. (a) Temperature gradient field values in log-10 scale for Arctic pack ice microstructure, with increasing connectivity from left to right (images courtesy of Don Perovich), with corresponding spectral functions μ(λ) displayed below in (b). Analogous to Fig. 3, as the fluid phase becomes increasingly connected on macroscopic length scales, a buildup of spectral measure mass at λ=0 shown in (b) leads to the formation of a δ function at λ=0; here, ϕ is the open-ocean area fraction for the images in (a). We see a corresponding switching in the values of the $m_i^{\mathrm{2}}$ from numerically zero, with the $m_i^{\mathrm{2}}\lesssim {\mathrm{10}}^{-\mathrm{30}}$ for the left and middle figure panels, to $m_i^{\mathrm{2}}\sim \mathrm{1}$ for the rightmost panel. The δ function at λ=1 is also analogous to that in Fig. 3. Like E₀ in Fig. 3, we take the average thermal gradient to be vertically oriented. Thermal conductivity is taken to be κ₁=0.57 W m⁻¹ K⁻¹ for ocean and κ₂=2.11 W m⁻¹ K⁻¹ for ice floes (Pringle et al., 2006).

From Theorem 1, the integral and χ₁E in Eqs. (5) and (6) have explicit representations in terms of the eigenvalues λ_i and eigenvectors w_i of M (Murphy et al., 2015, 2017a),

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& {\mathit{\chi }}_{\mathrm{1}}\mathit{E}=s\sum _i{\displaystyle \frac{m_i}{s-{\mathit{\lambda }}_i}}{\mathit{w}}_i,F\left(s\right)=\sum _i〈{\displaystyle \frac{m_i^{\mathrm{2}}}{s-{\mathit{\lambda }}_i}}〉,\\ & m_i={\mathit{\chi }}_{\mathrm{1}}{\mathit{w}}_i\cdot {\widehat{\mathit{e}}}_k,\end{array}\end{array}$$

where ${\widehat{\mathit{e}}}_k$ plays the role of a standard basis vector on the lattice. To compute μ, a non-standard generalization of the spectral theorem for matrices is required, due to the projective nature of the matrices χ₁ and Γ (Murphy et al., 2015). We developed a projection method that shows the spectral measure μ in Eq. (8) depends only on the eigenvalues and eigenvectors of random submatrices of Γ of size N₁≈ϕN. They correspond to diagonal components [χ₁]_ii=1, as the spectral weights m_i (Christoffel numbers) associated with eigenvectors satisfying χ₁w_i=0 are themselves zero, m_i=0 (Murphy et al., 2015). Fortunately, since these submatrices are much smaller for low volume fractions, this method greatly improves the efficiency and accuracy of numerical computations of μ.

The measure μ exhibits fascinating transitional behavior as a function of system connectivity. For example, in the case of a random resistor network (RRN) with a low volume fraction p of open bonds, as shown in Fig. 2a, there are spectrum-free regions at the spectral endpoints $\mathit{\lambda }=\mathrm{0},\mathrm{1}$ (Murphy and Golden, 2012; Murphy et al., 2015). However, as p approaches the percolation threshold p_c (Stauffer and Aharony, 1992; Torquato, 2002) and the system becomes increasingly connected, these spectral gaps shrink and then vanish (Murphy and Golden, 2012; Jonckheere and Luck, 1998), leading to the formation of δ components of μ at the spectral endpoints, precisely (Murphy and Golden, 2012) when p=p_c (and $p=\mathrm{1}-p_{\mathrm{c}}$ in d=3). This leads to critical behavior of σ^* for insulating/conducting and conducting/superconducting systems (Murphy and Golden, 2012). This gap behavior of μ has led (Golden, 1997c; Murphy and Golden, 2012) to a detailed description of these critical transitions in ${\mathit{\sigma }}^{*},$ which is analogous to the Lee–Yang–Ruelle–Baker description (Baker, 1990; Golden, 1997c) of the Ising model phase transition in the magnetization M. Moreover, using this gap behavior, all of the classical critical exponent scaling relations were recovered (Murphy and Golden, 2012; Golden, 1997c) without heuristic scaling forms (Efros and Shklovskii, 1976) but instead by using the rigorous integral representation for σ^* involving μ.

This spectral behavior emerges in all the systems mentioned above, such as the brine microstructure of sea ice (Golden et al., 1998a, 2007; Golden, 2009) as shown in Fig. 3, melt ponds on the surface of Arctic sea ice (Hohenegger et al., 2012) as shown in Fig. 4, and the sea ice pack itself (Murphy et al., 2017a) in Fig. 5. This also gives rise to critical behavior of the electric field, as shown in Fig. 3 for 2D cross sections of 3D brine microstructure, with E₀ 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 in the electric field's strength, with “hot spots” forming in geometric bottlenecks. A similar behavior is exhibited by the temperature gradient ∇T associated with the Stieltjes integral for the effective thermal conductivity κ^*, as shown for melt ponds atop Arctic sea ice in Fig. 4 and for Arctic pack ice in Fig. 5.

3.2 Generalization to rank deficient setting

In the periodic setting, for example, the matrix Laplacian is singular, so the matrix representation of $(-\mathrm{\Delta }{)}^{-\mathrm{1}}$ in Γ is not defined. We now extend the mathematical framework developed in Murphy et al. (2015) to this setting. To make the connection to the abstract Hilbert space (Golden and Papanicolaou, 1983) and full-rank matrix (Murphy et al., 2015) settings, we first give relevant details for these cases. Equation (6) for the abstract Hilbert space setting follows by applying the operator $-\mathbf{\nabla }(-\mathrm{\Delta }{)}^{-\mathrm{1}}$ to the formula $\mathbf{\nabla }\cdot \mathit{D}=\mathrm{0}$, yielding ΓD=0. Equation (6) then follows by using ΓE_f=E_f and ΓE₀=0 (Murphy et al., 2015), since E_f is in the range of Γ and E₀ is constant (Murphy et al., 2020a, 2017b, 2015). The matrix form of $\mathbf{\nabla }\cdot \mathit{D}=\mathrm{0}$ is $-{\mathbf{\nabla }}^T\mathit{D}=\mathrm{0}$, where ∇ now represents the finite-difference matrix representation of the gradient operator and −∇^T is the finite-difference representation of the divergence operator, with negative matrix Laplacian given by ∇^T∇ (Murphy et al., 2015). As before, in Murphy et al. (2015) we applied the matrix ∇(∇^T∇)⁻¹ to the formula $-{\mathbf{\nabla }}^T\mathit{D}=\mathrm{0}$, yielding ΓD=0, where $\mathbf{\Gamma }=\mathbf{\nabla }({\mathbf{\nabla }}^T\mathbf{\nabla }{)}^{-\mathrm{1}}{\mathbf{\nabla }}^T$, and Eq. (6) follows the same way as before, using the fact that E₀ is in the orthogonal complement of the range of ∇, so ΓE₀=0.

Now consider the singular value decomposition of the matrix gradient (Murphy et al., 2020a) of size m×n, say, ∇=UΣV^T. Here U is a m×n matrix satisfying U^TU=I_n; Σ is a n×n diagonal matrix with diagonal entries consisting of the singular values of ∇; and V is a n×n orthogonal matrix satisfying ${\mathbf{V}}^T\mathbf{V}={\mathbf{VV}}^T={\mathbf{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 Σ is an invertible matrix, and the matrix representation of Γ is given by Γ=UU^T. On the other hand, when the matrix gradient is singular we have $\mathbf{\Sigma }=\text{diag}({\mathbf{\Sigma }}_{\mathrm{1}},\mathrm{0},\mathrm{\dots },\mathrm{0})$, where the diagonal matrix Σ₁ contains the n₁ strictly positive singular values of Σ and the rest of the singular values have value 0. Denoting U₁ and V₁ to be the columns of U and V corresponding to the diagonal entries of Σ₁, we have $\mathbf{\nabla }={\mathbf{U}}_{\mathrm{1}}{\mathbf{\Sigma }}_{\mathrm{1}}{\mathbf{V}}_{\mathrm{1}}^T$, where Σ₁ is invertible and ${\mathbf{U}}_{\mathrm{1}}^T{\mathbf{U}}_{\mathrm{1}}={\mathbf{V}}_{\mathrm{1}}^T{\mathbf{V}}_{\mathrm{1}}={\mathbf{I}}_{n_{\mathrm{1}}}$. This enables us to write $-{\mathbf{\nabla }}^T\mathit{D}=\mathrm{0}$ as $-{\mathbf{V}}_{\mathrm{1}}{\mathbf{\Sigma }}_{\mathrm{1}}{\mathbf{U}}_{\mathrm{1}}^T\mathit{D}=\mathrm{0}$; hence ${\mathbf{U}}_{\mathrm{1}}^T\mathit{D}=\mathrm{0}$ and ${\mathbf{U}}_{\mathrm{1}}{\mathbf{U}}_{\mathrm{1}}^T\mathit{D}=\mathrm{0}$. Noting that the columns of U₁ span the range of the matrix gradient ∇, the matrix ${\mathbf{U}}_{\mathrm{1}}{\mathbf{U}}_{\mathrm{1}}^T$ is a projection onto the range of ∇ (Murphy et al., 2020a). Defining $\mathbf{\Gamma }={\mathbf{U}}_{\mathrm{1}}{\mathbf{U}}_{\mathrm{1}}^T$, Eq. (6) follows the same way as before, using the fact that E₀ is in the orthogonal complement of the range of ∇ so ΓE₀=0. This clearly generalizes the full-rank setting. More details will be published elsewhere.

4 Analytic continuation for polycrystalline media

Sea ice is a composite material with polycrystalline 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. 6a.

Our analysis of the transport properties of random, uniaxial polycrystalline media (Barabash and Stroud, 1999) in Gully et al. (2015), and a somewhat new formulation presented below, shows that the underlying mathematical framework is a direct analogue of that for two-phase random media discussed in Sect. 3. For simplicity, we discuss electrical permittivity ϵ, keeping in mind the broader applicability to thermal conductivity κ and electric conductivity σ, as well as viscoelastic modulus (Cherkaev, 2019), 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 Eq. (1), the local permittivity matrix of such media is given by (Milton, 2002; Barabash and Stroud, 1999):

$$\begin{array}{}\text{(10)}& {\displaystyle }\mathit{ϵ}(\mathit{x},\mathit{\omega })={\mathbf{R}}^T\mathbf{\Phi }\mathbf{R},\mathbf{\Phi }=\text{diag}({\mathit{ϵ}}_{\mathrm{1}},\mathrm{\dots },{\mathit{ϵ}}_d),\end{array}$$

where R(x,ω) is a random rotation matrix satisfying ${\mathbf{R}}^T={\mathbf{R}}^{-\mathrm{1}}$. For example, for d=2 we have

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathit{ϵ}={\mathbf{R}}^T\left[\begin{array}{cc}{\mathit{ϵ}}_{\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& {\mathit{ϵ}}_{\mathrm{2}}\end{array}\right]\mathbf{R},\mathbf{R}=\left[\begin{array}{rr}\cos \mathit{\theta }& -\sin \mathit{\theta }\\ \sin \mathit{\theta }& \cos \mathit{\theta }\end{array}\right],\end{array}$$

where $\mathit{\theta }=\mathit{\theta }(\mathit{x},\mathit{\omega })$ is the orientation angle, measured from the direction e₁ of the crystallite, which has an interior containing x∈ℝ^d for ω∈Ω. In higher dimensions, d≥3, the rotation matrix R is a composition of “basic” rotation matrices R_j, e.g., $\mathbf{R}={\prod }_{j=\mathrm{1}}^d{\mathbf{R}}_j$, where the matrix R_j(x,ω) rotates vectors in ℝ^d by an angle ${\mathit{\theta }}_j={\mathit{\theta }}_j(\mathit{x},\mathit{\omega })$ about the e_j axis. For example, in three dimensions

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& {\mathbf{R}}_{\mathrm{1}}=\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \cos {\mathit{\theta }}_{\mathrm{1}}& -\sin {\mathit{\theta }}_{\mathrm{1}}\\ \mathrm{0}& \sin {\mathit{\theta }}_{\mathrm{1}}& \cos {\mathit{\theta }}_{\mathrm{1}}\end{array}\right],\\ & {\mathbf{R}}_{\mathrm{2}}=\left[\begin{array}{ccc}\cos {\mathit{\theta }}_{\mathrm{2}}& \mathrm{0}& \sin {\mathit{\theta }}_{\mathrm{2}}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ -\sin {\mathit{\theta }}_{\mathrm{2}}& \mathrm{0}& \cos {\mathit{\theta }}_{\mathrm{2}}\end{array}\right],\\ & {\mathbf{R}}_{\mathrm{3}}=\left[\begin{array}{ccc}\cos {\mathit{\theta }}_{\mathrm{3}}& -\sin {\mathit{\theta }}_{\mathrm{3}}& \mathrm{0}\\ \sin {\mathit{\theta }}_{\mathrm{3}}& \cos {\mathit{\theta }}_{\mathrm{3}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right].\end{array}\end{array}$$

In the case of uniaxial polycrystalline media, the local permittivity along one of the crystal axes has the value ϵ₁, while the permittivity along all the other crystal axes has the value ϵ₂, so $\mathbf{\Phi }=\text{diag}({\mathit{ϵ}}_{\mathrm{1}},{\mathit{ϵ}}_{\mathrm{2}})$ for 2D (which is the general setting for 2D) and $\mathbf{\Phi }=\text{diag}({\mathit{ϵ}}_{\mathrm{1}},{\mathit{ϵ}}_{\mathrm{2}},{\mathit{ϵ}}_{\mathrm{2}})$ for 3D. Equation (10) can be written in a more suggestive form in terms of the matrix $\mathbf{C}=\text{diag}(\mathrm{1},\mathrm{0},\mathrm{\dots },\mathrm{0})$:

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathit{ϵ}(\mathit{x},\mathit{\omega })={\mathit{ϵ}}_{\mathrm{1}}{\mathbf{X}}_{\mathrm{1}}(\mathit{x},\mathit{\omega })+{\mathit{ϵ}}_{\mathrm{2}}{\mathbf{X}}_{\mathrm{2}}(\mathit{x},\mathit{\omega }),\end{array}$$

which is an analogue of Eq. (2). Here X₁=R^TCR and ${\mathbf{X}}_{\mathrm{2}}={\mathbf{R}}^T(\mathbf{I}-\mathbf{C})\mathbf{R}$, where I is the identity matrix on ℝ^d. Since ${\mathbf{R}}^T={\mathbf{R}}^{-\mathrm{1}}$ and C is a diagonal projection matrix satisfying C²=C, it is clear that the X_i, $i=\mathrm{1},\mathrm{2}$, are mutually orthogonal projection matrices satisfying

$$\begin{array}{}\text{(14)}& {\displaystyle }{\mathbf{X}}_j^T={\mathbf{X}}_j,{\mathbf{X}}_j{\mathbf{X}}_k={\mathbf{X}}_j{\mathit{\delta }}_{jk},{\mathbf{X}}_{\mathrm{1}}+{\mathbf{X}}_{\mathrm{2}}=\mathbf{I},\end{array}$$

which are also properties of the characteristic functions χ_j in Sect. 3.

Figure 6 Spectral analysis of polycrystalline media. (a) Cross sections of polycrystalline microstructure for granular and columnar sea ice (Jean-Louis Tyson). (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 μ(λ), shown along with its theoretical prediction for such isotropic media (Milton, 2002). (d) An example value of the complex effective permittivity of isotropic polycrystalline media captured by first- and second-order bounds (Gully et al., 2015).

Equations (3) and (4) are also satisfied in this polycrystalline setting (Golden and Papanicolaou, 1983). Similar to the derivation of Eq. (6) in Sect. 3, a resolvent representation for X₁E follows by applying the operator $-\mathbf{\nabla }(-\mathrm{\Delta }{)}^{-\mathrm{1}}$ to the formula $\mathbf{\nabla }\cdot \mathit{D}=\mathrm{0}$, yielding ΓD=0. Then, using ΓE_f=E_f and ΓE₀=0 (Murphy et al., 2015) yields the following analogue of Eq. (6):

$$\begin{array}{}\text{(15)}& {\displaystyle }{\mathbf{X}}_{\mathrm{1}}\mathit{E}=s(sI-G{)}^{-\mathrm{1}}{\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_k,G={\mathbf{X}}_{\mathrm{1}}\mathrm{\Gamma }{\mathbf{X}}_{\mathrm{1}},\end{array}$$

yielding the integral representation in Eq. (5) for $F\left(s\right)=\langle \left[\right(sI-G{)}^{-\mathrm{1}}{\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_k]\cdot {\mathit{e}}_k\rangle$. As in the two-component setting, a critical feature of Eq. (5) is that the component parameters in s are separated from the geometrical information in μ. Information about the geometry enters through the moments in Eq. (7), with G given in Eq. (15) and χ₁ replaced by X₁. The mass of the measure μ_jk is given by

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathit{\mu }}_{jk}^{\mathrm{0}}=\langle {\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_j\cdot {\mathit{e}}_k\rangle ,{\mathit{\mu }}_{kk}^{\mathrm{0}}=\langle |{\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_k{|}^{\mathrm{2}}\rangle ,\end{array}$$

where the second equality follows from the fact that X₁ is a real-symmetric projection matrix. The statistical average $\langle |{\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_k{|}^{\mathrm{2}}\rangle$ in Eq. (16) can be thought of as the “mean orientation” or as the percentage of crystallites oriented in the kth direction. For example, in the case of two-dimensional polycrystalline media, d=2, Eq. (11) implies that

$$\begin{array}{}\text{(17)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{11}}^{\mathrm{0}}=\langle {\cos }^{\mathrm{2}}\mathit{\theta }\rangle ,{\mathit{\mu }}_{\mathrm{22}}^{\mathrm{0}}=\langle {\sin }^{\mathrm{2}}\mathit{\theta }\rangle ,{\mathit{\mu }}_{\mathrm{12}}^{\mathrm{0}}=\langle \sin \mathit{\theta }\cos \mathit{\theta }\rangle .\end{array}$$

Generalizing Eq. (12), with $\mathbf{R}={\prod }_{j=\mathrm{1}}^d{\mathbf{R}}_j$, to dimensions d≥3 shows that ${\mathit{\mu }}_{jk}^{\mathrm{0}}$ is a linear combination of averages of the form $\langle {\prod }_i{\cos }^{n_i}{\mathit{\theta }}_i{\sin }^{m_i}{\mathit{\theta }}_i\rangle$, where $n_i,m_i=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots }$.

Given partial information on the microgeometry incorporated into the moments μ_n, the integral representation (5) for this polycrystalline setting yields rigorous forward bounds for the effective parameters of composites (Gully et al., 2015; Milton, 2002), as shown in Fig. 6d. The integral representation can also be used to obtain 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 (Gully et al., 2015; Milton, 2002); see Sect. 5 for more details.

Spectral measure computations for uniaxial polycrystalline materials

Computing the spectral measure μ for a given polycrystalline microgeometry first involves discretizing the composite into a square lattice with vertex values in the range [0,2π] corresponding to the crystallite orientation angles at each vertex location. On this square lattice the action of the differential operators ∇ and ∇⋅ are defined in terms of forward- and backward-difference operators (Golden, 1992; Murphy et al., 2015). Then the key operator X₁ΓX₁, which depends on the geometry of the network via X₁, becomes a real-symmetric matrix M. Here Γ is as in Sect. 3.1, and X₁ is a banded (random) projection matrix that determines the geometry of the polycrystalline medium. In this setting, the integral and X₁E in Eqs. (5) and (6) have explicit representations in terms of the eigenvalues λ_i and eigenvectors w_i of M (Murphy et al., 2015) given by Eq. (9), and similarly the spectral measure is given by Eq. (8), with χ₁ replaced by X₁.

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 1 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.

Theorem 2.

For each ω∈Ω, let M(ω)=W(ω)Λ(ω) W^T(ω) be the eigenvalue decomposition of the real-symmetric matrix M(ω)=X₁(ω) Γ X₁(ω). Here, the columns of the matrix W(ω) consist of the orthonormal eigenvectors w_i(ω), $i\mathrm{=}\mathrm{1}\mathrm{,}\mathrm{\dots }\mathrm{,}N$, of M(ω), and the diagonal matrix $\mathbf{\Lambda }\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\text{diag}\mathrm{\left(}{\mathit{\lambda }}_{\mathrm{1}}\mathrm{\right(}\mathit{\omega }\mathrm{)}\mathrm{,}\mathrm{\dots }\mathrm{,}{\mathit{\lambda }}_N\mathrm{(}\mathit{\omega }\mathrm{\left)}\mathrm{\right)}$ involves its eigenvalues λ_i(ω). Denote ${\mathbf{Q}}_i\mathrm{=}{\mathit{w}}_i{\mathit{w}}_i^T$ the projection matrix onto the eigenspace spanned by w_i. The electric field E(ω) satisfies $\mathit{E}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}{\mathit{E}}_{\mathrm{0}}\mathrm{+}{\mathit{E}}_f\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$, with E₀=〈E(ω)〉 and ΓE(ω)=E_f(ω), and the effective complex permittivity tensor ϵ^* has components ${\mathit{ϵ}}_{jk}^{\mathrm{*}}$, $j\mathrm{,}k\mathrm{=}\mathrm{1}\mathrm{,}\mathrm{\dots }\mathrm{,}d$, which satisfy

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& {\mathit{ϵ}}_{jk}^{*}={\mathit{ϵ}}_{\mathrm{2}}({\mathit{\delta }}_{jk}-F_{jk}(s\left)\right),F_{jk}\left(s\right)=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{\mathrm{d}{\mathit{\mu }}_{jk}\left(\mathit{\lambda }\right)}{s-\mathit{\lambda }}},\\ & \mathrm{d}{\mathit{\mu }}_{jk}\left(\mathit{\lambda }\right)=\sum _{i=\mathrm{1}}^N\langle {\mathit{\delta }}_{{\mathit{\lambda }}_i}\left(\mathrm{d}\mathit{\lambda }\right){\mathbf{X}}_{\mathrm{1}}{\mathbf{Q}}_i{\widehat{\mathit{e}}}_j\cdot {\widehat{\mathit{e}}}_k\rangle .\end{array}\end{array}$$

To numerically compute μ, a non-standard generalization of the spectral theorem for matrices is required due to the projective nature of the matrices X₁ and Γ (Murphy et al., 2015). In particular, a projection method analogous to that in Murphy et al. (2015) shows that the spectral measure μ in Eq. (18) depends only on the eigenvalues and eigenvectors of the upper left N₁×N₁ block of the matrix RΓR^T, where $N_{\mathrm{1}}=N/d$. These submatrices are smaller by a factor of d, which decreases the numerical cost of computations of μ by a factor of d³. In Fig. 6, computations of the displacement field D 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 ϵ^* is diagonal, such as the setting of isotropically oriented crystallites, the spectral measure for an infinite system is known in closed form (Milton, 2002) to be $\mathrm{d}\mathit{\mu }\left(\mathit{\lambda }\right)=\left(\sqrt{(\mathrm{1}-\mathit{\lambda })/\mathit{\lambda }}\right)(\mathrm{d}\mathit{\lambda }/\mathit{\pi })$, as shown in Fig. 6c. This measure has a singularity at λ=0, which indicates that the material is electrically conductive on macroscopic length scales (Murphy et al., 2015; Murphy and Golden, 2012). When the polycrystalline material is isotropic, both the mass and first moment of the measure μ are known, which enables two nested bounds for ϵ to be computed (Gully et al., 2015), as shown in Fig. 6d.

5 Inverse homogenization: recovery of information about the microstructure of composites

Developed originally for the effective complex permittivity ϵ^*, the integral representation (Eq. 5) yields rigorous forward bounds for the effective permittivity ϵ^* of two-component composites formed of materials with permittivity ϵ₁ and ϵ₂, given partial information on the microgeometry via the moments μ_n (Bergman, 1980; Milton, 1980; Bergman, 1982; Golden and Papanicolaou, 1983). One can also use the integral representation to recover information about the structure of the composite material – this is the problem of inverse homogenization (Cherkaev, 2001). For inverse homogenization, it is important that the representation (Eq. 5) separates information about the properties of the phases contained in the parameter s from information about the microgeometry contained in the measure μ and its moments ${\mathit{\mu }}_n=\langle G^n{\mathit{\chi }}_{\mathrm{1}}{\mathit{e}}_k\cdot {\mathit{\chi }}_{\mathrm{1}}{\mathit{e}}_k\rangle$ in Eq. (7) via higher-order correlation functions of the geometry function χ₁.

In principle, the spectral measure μ and its moments μ_n contain all the geometrical information about the composite. For example, the mass μ₀ is the volume fraction ϕ of the first component in the composite,

$$\begin{array}{}\text{(19)}& {\mathit{\mu }}_{\mathrm{0}}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathrm{d}\mathit{\mu }\left(z\right)=\langle {\mathit{\chi }}_{\mathrm{1}}\rangle =\mathit{\varphi },\end{array}$$

and the fraction of the second phase is 1−ϕ. Connectivity information is also embedded in the spectral measure.

The basis for inverse homogenization is provided by the uniqueness theorem (Cherkaev, 2001), which formulates the conditions under which the measure μ in the representation (Eq. 5) can be uniquely reconstructed from measured data. For instance, complex permittivity data measured for a range of frequencies of the applied electromagnetic field are sufficient to uniquely recover the measure μ in Eq. (5) (Cherkaev, 2001). Such data are also sufficient for the unique and stable reconstruction of the moments μ_n (Cherkaev and Ou, 2008), provided the permittivity of one of the phases is frequency dependent. Two major approaches to inverse homogenization are (i) reconstruction of the measure µ (Cherkaev, 2001; Cherkaev and Ou, 2008; Day and Thorpe, 1996; Zhang and Cherkaev, 2009; Bonifasi-Lista and Cherkaev, 2009; Bonifasi-Lista et al., 2009; Cherkaev and Bonifasi-Lista, 2011; Day and Thorpe, 1999; Day et al., 2000; Golden et al., 2011; Cherkaev, 2020) and then calculating its moments and (ii) inverse bounds for the structural parameters, such as the volume fractions of the components (McPhedran et al., 1982; McPhedran and Milton, 1990; Cherkaev and Tripp, 1996; Cherkaev and Golden, 1998; Cherkaev, 2001; Cherkaev and Ou, 2008), orientation of the crystals (Gully et al., 2015), or connectedness (Orum et al., 2012) of one phase.

When only a few data points are available, though the uniqueness theorem (Cherkaev, 2001) is not immediately applicable, one can outline a set of measures consistent with the measurements,

$$\begin{array}{}\text{(20)}& \mathcal{M}=\mathit{\{}\mathit{\mu }:F_{\mathit{\mu }}(s)=\mathrm{1}-{\mathit{ϵ}}^{\ast }/{\mathit{ϵ}}_{\mathrm{2}}\mathit{\}},\end{array}$$

and then determine an interval confining the first moment of the measure μ 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 (Cherkaev and Tripp, 1996; Cherkaev and Golden, 1998; Tripp et al., 1998). When the requirements for the measurements needed to uniquely reconstruct the spectral measure μ established by the uniqueness theorem are satisfied, the set ℳ 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 ${\mathit{\mu }}_{\mathrm{0}},{\mathit{\mu }}_{\mathrm{1}},\mathrm{\dots }$ are the S-equivalent structures (Cherkaev, 2001) that are not distinguishable by homogenized measurements.

An equivalent representation for the function F(s) in Eq. (5) using a logarithmic potential of the measure μ in the complex s plane is (Cherkaev, 2001)

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& F\left(s\right)={\displaystyle \frac{\mathrm{d}}{\mathrm{d}s}}\int \ln |s-z|\mathrm{d}\mathit{\mu }\left(z\right),\\ & \mathrm{d}/\mathrm{d}s=(\partial /\partial x-i\partial /\partial y),s={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{ϵ}}_{\mathrm{1}}/{\mathit{ϵ}}_{\mathrm{2}}}}.\end{array}\end{array}$$

The solution to the inverse problem of recovering the measure μ is constructed by solving the minimization problem:

$$\begin{array}{}\text{(22)}& {\text{min}}_{\mathit{\mu }}\left|\right|A\mathit{\mu }-F|{|}^{\mathrm{2}},F(s)=\mathrm{1}-{\mathit{ϵ}}^{*}(s)/{\mathit{ϵ}}_{\mathrm{2}}.\end{array}$$

Here A is the integral operator in Eq. (21) or in Eq. (5); the norm is the L² norm; F is the given function of the measured data; and s∈𝒞, where 𝒞 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⁻¹ leads to arbitrarily large variations in the solution, and the problem requires regularization to design a stable numerical algorithm (Cherkaev, 2001). Regularized inversion schemes and stable reconstruction algorithms to recover μ and its moments from data on the effective complex permittivity were developed in Cherkaev (2001, 2004), Cherkaev and Ou (2008), Bonifasi-Lista and Cherkaev (2009), and Cherkaev and Bonifasi-Lista (2011). They are based on L², total variation (TV), non-negativity constraints, and constrained Padé approximations of the measure μ (Zhang and Cherkaev, 2009). In applications to imaging of bone structure, spectral measures μ are computed with regularization algorithms based on L² constrained minimization, from electromagnetic (Bonifasi-Lista and Cherkaev, 2009; Cherkaev and Bonifasi-Lista, 2011; Golden et al., 2011) and viscoelastic (Bonifasi-Lista and Cherkaev, 2008; Bonifasi-Lista et al., 2009; Cherkaev and Bonifasi-Lista, 2011) 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.

The first application of Stieltjes representations to the elastic properties of two-phase composites can be found in Kantor and Bergman (1982, 1984). Using hydrostatic and deviatoric projections Λ_h and Λ_s onto the orthogonal subspaces of the second-order tensors which are proportional to the identity tensor and trace-free, the Stieltjes integral representation was generalized in Cherkaev and Bonifasi-Lista (2011) to the effective viscoelastic modulus and to two-dimensional viscoelastic polycrystalline materials (Cherkaev, 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 (Bonifasi-Lista and Cherkaev, 2008; Cherkaev and Bonifasi-Lista, 2011; Cherkaev, 2020) for successful recovery of the volume fractions of the phases in a composite from a known effective viscoelastic shear modulus.

Other approaches to the volume fraction bounds include Engström (2005), Milton (2012), and Thaler and Milton (2014), based on estimates for higher-order moments and on variational bounds, as well as direct inversion of known formulas or mixing rules (Bergman and Stroud, 1992; Levy and Cherkaev, 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.

Spectral coupling of various properties of composites

An important application of inverse homogenization is for indirectly evaluating material properties through cross-coupling (Cherkaev, 2001; Milton, 2002). Different properties of composites are coupled through their microgeometry; this relationship can be used for estimating properties which are difficult to measure from other effective properties. The conventional approaches are based on empirical and semi-empirical relations, such as, for instance, Kozeny–Carman or Katz–Thompson. These relations estimate the permeability of a porous material characterizing the microstructure by a ”formation factor” F, which relates the properties of one phase in the composite to the effective properties of the material (Sahimi, 1995; Torquato, 2002; Wong et al., 1984; Wong, 1988).

In the spectral coupling method (Cherkaev, 2001) based on the Stieltjes representation (5), the spectral measure μ associated with the geometric structural function, couples various properties of the same material. Spectral coupling (Cherkaev, 2001, 2004; Cherkaev and Zhang, 2003; Cherkaev and Bonifasi-Lista, 2011) 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 (Cherkaev and 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 Cherkaev and Bonifasi-Lista (2011), 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.

5.1 Bounds for the moments of the spectral measure

The second approach to the inverse homogenization problem is calculating inverse bounds for the structural parameters, such as the volume fraction of each of the components (McPhedran et al., 1982; McPhedran and Milton, 1990; Cherkaev and Tripp, 1996; Cherkaev and Golden, 1998; Cherkaev, 2001), orientation of the crystals (Gully et al., 2015), or connectedness (Orum et al., 2012) of, say, the conducting phase. An analytical approach to estimating the volume fractions of materials in a composite (Cherkaev and Tripp, 1996; Cherkaev and Golden, 1998; Tripp et al., 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 (Cherkaev and Golden, 1998).

The inverse bounds are derived using analyticity of the effective complex permittivity of the composite. The first-order bounds $p_l^{\left(\mathrm{1}\right)}$ and $p_u^{\left(\mathrm{1}\right)}$ for the volume fraction ϕ give the lower and upper bounds for the zeroth moment μ₀ of the measure μ or its mass in Eq. (19) (Cherkaev and Tripp, 1996; Cherkaev and Golden, 1998):

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& p_l^{\left(\mathrm{1}\right)}\le \mathit{\varphi }\le p_u^{\left(\mathrm{1}\right)},p_l^{\left(\mathrm{1}\right)}=|f{|}^{\mathrm{2}}{\displaystyle \frac{\text{Im}\left(\overline{s}\right)}{\text{Im}\left(f\right)}},\\ & p_u^{\left(\mathrm{1}\right)}=\mathrm{1}-{\displaystyle \frac{\left|g{|}^{\mathrm{2}}\text{Im}\right(\overline{t})}{\text{Im}\left(g\right)}}.\end{array}\end{array}$$

Here, $t=\mathrm{1}-s$, f is the known value of F(s), and g is the known value of $G\left(t\right)=\mathrm{1}-{\mathit{ϵ}}^{*}/{\mathit{ϵ}}_{\mathrm{1}}$.

Figure 7 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 (Cherkaev and Golden, 1998). The small lens-shaped domains each contain ϵ^* 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 (Orum et al., 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 (Gully et al., 2015) for the permittivity of sea ice (left) together with the measured effective permittivity of sea ice in (Arcone et al., 1986). 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.

First- and second-order forward and inverse bounds are illustrated in Fig. 7a (Cherkaev and Golden, 1998), where first-order bounds for the effective complex permittivity of all anisotropic composites that could be formed from two materials of permittivity ϵ₁ and ϵ₂ 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 ϕ of the first component equal to $p_l^{\left(q\right)}$ and $p_u^{\left(q\right)}$, $q=\mathrm{1},\mathrm{2}$. The points $p_l^{\left(q\right)}$ and $p_u^{\left(q\right)}$ 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 ϵ^*(k), $k=\mathrm{1},\mathrm{\dots },N$ corresponding to composites with the same structure, the bounds for the fraction ϕ of the first phase in the material are given by an intersection of all admissible intervals $p_l^{\left(q\right)}\left(k\right)\le \mathit{\varphi }\le p_u^{\left(q\right)}\left(k\right)$ (Cherkaev and Tripp, 1996):

$$\begin{array}{}\text{(24)}& \begin{array}{rl}& P_l^{\left(q\right)}=\underset{k}{max}{p}_l^{\left(q\right)}\left(k\right)\le \mathit{\varphi }\le \underset{k}{min}{p}_u^{\left(q\right)}\left(k\right)=P_u^{\left(q\right)},\\ & q=\mathrm{1},\mathrm{2}.\end{array}\end{array}$$

Here, $p_l^{\left(q\right)}\left(k\right)$ and $p_u^{\left(q\right)}\left(k\right)$ are lower and upper bounds for the volume fraction ϕ calculated using the effective complex permittivity ϵ^∗(k), respectively, and q is the order of the bounds, with q=1 for a general mixture and q=2 for an isotropic composite.

In Cherkaev and Golden (1998), 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 ϵ^* of sea ice (Arcone et al., 1986) at −6 and −11 ^∘C with fractions of brine ϕ=0.036 and ϕ=0.0205. Sea ice was considered 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 $\mathrm{0.0213}\le \mathit{\varphi }\le \mathrm{0.0664}$ and $\mathrm{0.0119}\le \mathit{\varphi }\le \mathrm{0.0320}$, for the data set 1 and 2, respectively. The second-order inverse bounds derived with the assumption of 2D isotropy in the horizontal plane give the following estimates for the brine volume fraction: $\mathrm{0.0333}\le \mathit{\varphi }\le \mathrm{0.0422}$ for the first data set with brine volume ϕ=0.036 and $\mathrm{0.0189}\le \mathit{\varphi }\le \mathrm{0.0213}$ for the second data set with volume fraction of brine ϕ=0.0205. The first-order bounds are further extended to polycrystalline materials and allow estimation of the mean crystal orientation (Gully et al., 2015).

5.2 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 (Bruno, 1991). For a matrix-particle composite with separated inclusions, tighter bounds on the effective complex permittivity may be obtained. In Orum et al. (2012), sea ice is considered 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_{\mathrm{b}}/r_{\mathrm{i}}$. The minimal separation of brine inclusions is $\mathrm{2}(r_{\mathrm{i}}-r_{\mathrm{b}})=\mathrm{2}{r}_{\mathrm{i}}(\mathrm{1}-q)$. In this case, as it is shown in Bruno (1991), the support of μ in Eq. (5) lies in an interval [s_m,s_M], $\mathrm{0}<s_m<s_M<\mathrm{1}$ such that $s_m=\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{1}-q^{\mathrm{2}})$ and $s_M=\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{1}+q^{\mathrm{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. 7), 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 ϕ 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 Orum et al. (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. 7 (Orum et al., 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.

5.3 Extension to polycrystalline composites

The method of inverse bounds (Cherkaev and Tripp, 1996; Cherkaev and Golden, 1998; Tripp et al., 1998) for structural parameters of a composite from measured effective properties was extended to polycrystalline materials in Gully et al. (2015). In the case of a uniaxial polycrystalline composite, Gully et al. (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 (Weeks and Ackley, 1982), 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, which affects melt pond evolution, nutrient replenishment, brine convection, and other mesoscale processes in the ice cover.

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 Gully et al. (2015) are shown in Fig. 7c. 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 ${\mathit{\mu }}_{kk}^{\mathrm{0}}$ in Eq. (16) of the measure μ in the integral representation of the effective properties of a uniaxial polycrystalline material is ${\mathit{\mu }}_{kk}^{\mathrm{0}}=\langle |{\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_k{|}^{\mathrm{2}}\rangle$, where the statistical average $\langle |{\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_k{|}^{\mathrm{2}}\rangle$ can be viewed as the “mean crystal orientation” related to the percentage of crystallites oriented in the kth direction. Extending the inverse bounds method (Cherkaev and Tripp, 1996; Cherkaev and Golden, 1998; Tripp et al., 1998) to polycrystalline materials, the inverse polycrystalline bounds (Gully et al., 2015) estimate the mean crystal orientation by bounding the zeroth moment ${\mathit{\mu }}_{kk}^{\mathrm{0}}$ of the measure μ 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 Eq. (23):

$$\begin{array}{}\text{(25)}& \begin{array}{rl}& \langle {\mathit{e}}_k^T{\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_k{\rangle }_l\le \langle {\mathit{e}}_k^T{\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_k\rangle \le \langle {\mathit{e}}_k^T{\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_k{\rangle }_u,\\ & \langle {\mathit{e}}_k^T{\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_k{\rangle }_l=|f{|}^{\mathrm{2}}{\displaystyle \frac{\text{Im}\left(\stackrel{\mathrm{‾}}{s}\right)}{\text{Im}\left(f\right)}},\\ & \langle {\mathit{e}}_k^T{\mathbf{X}}_{\mathrm{1}}{\mathit{e}}_k{\rangle }_u=\mathrm{1}-|g{|}^{\mathrm{2}}{\displaystyle \frac{\text{Im}\left(\stackrel{\mathrm{‾}}{t}\right)}{\text{Im}\left(g\right)}},\end{array}\end{array}$$

Here, X₁ is defined in the polycrystalline section as X₁=R^TCR; f is the known value of F(s); and g is the known value of $G\left(t\right)=\mathrm{1}-{\mathit{ϵ}}^{*}/{\mathit{ϵ}}_{\mathrm{1}}$, with $t=\mathrm{1}-s$.

Inverse polycrystalline bounds computed in Gully et al. (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. 6a (top), the inverse mean crystal orientation bounds (Gully et al., 2015) estimate the deviation angle as $\mathit{\pi }/\mathrm{2}\pm .02$ (with the true value $\mathit{\pi }/\mathrm{2}$). The inverse mean crystal orientation bounds (Gully et al., 2015) for columnar ice (see Fig. 6a, bottom) estimate the angle of deviation of the crystal's axis from the vertical as 20±8^∘. 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.

6 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 (Washington and Parkinson, 1986) and Earth's ecosystems (Di Lorenzo et al., 2013). Advection of geophysical fluids intensifies the dispersion and large-scale transport of heat (Moffatt, 1983), pollutants (Csanady, 1963; Beychok, 1994; Samson, 1988), and nutrients (Di Lorenzo et al., 2013; Hofmann and Murphy, 2004) diffusing in their environment. In sea ice dynamics, where the ice cover couples the atmosphere to the polar oceans (Washington and Parkinson, 1986), the transport of sea ice can also be enhanced by eddy fluxes and large-scale coherent structures in the ocean (Watanabe and Hasumi, 2009; Lukovich et al., 2015). 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 (Thomas and Dieckmann, 2003; Golden et al., 2007; Golden, 2009). Convective brine flow through the porous microstructure can enhance thermal transport through the sea ice layer (Lytle and Ackley, 1996; Worster and Jones, 2015).

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 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 with an effective diffusion tensor, which depends on the geometry of the velocity field and the local diffusivity. In Taylor (1921), it was first shown that a long-time, large-scale dispersion of passive scalars can be described by an effective diffusivity tensor κ^*. Motivated by Papanicolaou and Varadhan (1982), 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 (McLaughlin et al., 1985). The effective parameter problem of (anomalous) super-diffusion and sub-diffusion is given in Biferale et al. (1995) and Fannjiang (2000). Based on McLaughlin et al. (1985), Avellaneda and Majda (1989, 1991) adapted the ACM (Golden and Papanicolaou, 1983) to the advection diffusion equation and obtained a Stieltjes integral representation of the effective diffusivity tensor κ^* for flows with zero mean drift involving the Péclet number ξ of the flow. This representation encapsulates the geometric complexity of the flow in a spectral measure associated with a random Hermitian operator (or matrix). Using methods developed for composite media (Milton, 2002), they obtained rigorous bounds on the components of κ^*. Moreover, in analogy to methods developed for composites (Milton, 2002), they also found velocity fields which realize these bounds, such as the famous confocal sphere configurations that realize the Hashin–Shtrikman bounds of composites (Hashin and Shtrikman, 1962; Avellaneda and Majda, 1991). Remarkably, this method has also been extended to time-dependent flows (Avellaneda and Vergassola, 1995; Murphy et al., 2017b); flows with incompressible nonzero effective drift (Pavliotis, 2002; Fannjiang and Papanicolaou, 1994); flows where particles diffuse according to linear collisions (Pavliotis, 2010); and solute transport in porous media (Bhattacharya, 1999), 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 κ^*.

We now briefly describe our recent results on this framework (Murphy et al., 2017b, 2020a). It is an important example of how Stieltjes integral representations can provide a rigorous basis for analysis of problems for sea ice involving advection diffusion processes. The dispersion of a cloud of passive scalars with density ϕ(t,x) (not to be confused with earlier usage as the brine volume fraction) diffusing with molecular diffusivity ε and being advected by an incompressible velocity field u(t,x) satisfying $\mathbf{\nabla }\cdot \mathit{u}=\mathrm{0}$ is described by the advection–diffusion equation as follows:

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}=\mathit{u}\cdot \mathbf{\nabla }\mathit{\varphi }+\mathit{\epsilon }\mathrm{\Delta }\mathit{\varphi },\mathit{\varphi }(\mathrm{0},\mathit{x})={\mathit{\varphi }}_{\mathrm{0}}\left(\mathit{x}\right).\end{array}$$

Here, the initial density ϕ₀(x) and the fluid velocity field u are assumed to be given. In Eq. (26), the molecular diffusion constant is ε>0, and $\mathrm{\Delta }=\mathbf{\nabla }\cdot \mathbf{\nabla }={\mathrm{\nabla }}^{\mathrm{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 ε. To simplify our presentation, we assume that the velocity field u in Eq. (26) is temporally and spatially periodic. Non-dimensionalization (Murphy et al., 2020a) and homogenization (McLaughlin et al., 1985) of Eq. (26) shows that long-time, large-volume (or area) macroscopic thermal transport is described by a diffusion equation involving an averaged scalar density $\overline{\mathit{\varphi }}$ and a symmetric, constant (Pavliotis, 2002) effective diffusivity tensor κ^* (McLaughlin et al., 1985),

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\partial \overline{\mathit{\varphi }}(t,\mathit{x})}{\partial t}}=\mathbf{\nabla }\cdot \left[{\mathit{\kappa }}^{*}\mathbf{\nabla }\overline{\mathit{\varphi }}\right(t,\mathit{x}\left)\right],\overline{\mathit{\varphi }}(\mathrm{0},\mathit{x})={\mathit{\varphi }}_{\mathrm{0}}\left(\mathit{x}\right).\end{array}$$

For simplicity, we consider a diagonal coefficient ${\mathit{\kappa }}_{kk}^{*}$, $k=\mathrm{1},\mathrm{\dots },d$, of κ^*; set ${\mathit{\kappa }}^{*}={\mathit{\kappa }}_{kk}^{*}$; and write u=u₀v, where u₀ is the (constant) strength of u and v is a non-dimensional velocity field containing the geometric and dynamic information about u. In these non-dimensional variables, the Péclet number ξ and molecular diffusivity ε are related by $\mathit{\xi }=\mathrm{1}/\mathit{\epsilon }$ (Murphy et al., 2017b, 2020a).

Using a mathematical framework that is strikingly similar to that in Sect. 3, the effective diffusivity has the following Stieltjes integral representation (McLaughlin et al., 1985; Avellaneda and Majda, 1991; Murphy et al., 2017b, 2020a):

$$\begin{array}{}\text{(28)}& {\displaystyle }{\mathit{\kappa }}^{*}=\mathit{\epsilon }(\mathrm{1}+\langle |\mathbf{\nabla }{\mathit{\psi }}_k{|}^{\mathrm{2}}\rangle ),\langle |\mathbf{\nabla }{\mathit{\psi }}_k{|}^{\mathrm{2}}\rangle =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}\mathit{\nu }\left(\mathit{\lambda }\right)}{{\mathit{\epsilon }}^{\mathrm{2}}+{\mathit{\lambda }}^{\mathrm{2}}}},\end{array}$$

where 〈⋅〉 denotes averaging over the space–time period cell for periodic flows (Murphy et al., 2017b, 2020a) or statistical average for random flows (Avellaneda and Majda, 1989; Avellaneda and Vergassola, 1995), and ψ_k is the solution to a cell problem (McLaughlin et al., 1985; Murphy et al., 2017b). An equivalent statement which emphasizes the connection to the two-component composites setting in Eq. (5) is

$$\begin{array}{}\text{(29)}& F\left(\mathit{\epsilon }\right)=\mathrm{1}-{\displaystyle \frac{{\mathit{\kappa }}^{*}}{\mathit{\epsilon }}}=-\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}\mathit{\nu }\left(\mathit{\lambda }\right)}{{\mathit{\epsilon }}^{\mathrm{2}}+{\mathit{\lambda }}^{\mathrm{2}}}}.\end{array}$$

The vector field $\mathit{E}(t,\mathit{x})=\mathbf{\nabla }{\mathit{\psi }}_k(t,\mathit{x})+{\mathit{e}}_k$ satisfies Eq. (3) for two-component composite materials, with D=ϵE, $\mathit{ϵ}=\mathit{\epsilon }\mathbf{I}+\mathbf{S}$, $\mathbf{S}=(-\mathbf{\Delta }{)}^{-\mathrm{1}}{\partial }_t+\mathbf{H}$, and ϵ plays the role of the medium's electrical permittivity tensor (Murphy et al., 2017b, 2020a). Here, ∂_t denotes partial differentiation with respect to time, and H(t,x) is the stream matrix, given in terms of the incompressible velocity field $\mathit{v}=\mathbf{\nabla }\cdot \mathbf{H}$, that satisfies ${\mathbf{H}}^T=-\mathbf{H}$ (Avellaneda and Majda, 1991, 1989). When the flow is time-independent, v=v(x), then ψ_k=ψ_k(x) and S=H(x). Moreover ${\mathit{\kappa }}^{*}={\mathit{ϵ}}^{*}$, with ${\mathit{ϵ}}^{*}=({\mathit{ϵ}}^{*}{)}_{kk}$ defined above (Murphy et al., 2017b). The integral representation for κ^* in Eq. (28) follows from the resolvent formula

$$\begin{array}{}\text{(30)}& {\displaystyle }\mathbf{\nabla }{\mathit{\psi }}_k=(\mathit{\epsilon }\mathbf{I}+\mathrm{\Gamma }\mathbf{S}\mathrm{\Gamma }{)}^{-\mathrm{1}}{\mathit{g}}_k,{\mathit{g}}_k=-\mathrm{\Gamma }\mathbf{H}{\mathit{e}}_k,\end{array}$$

which is an analogue of Eq. (6). The operator ΓSΓ is antisymmetric due to the asymmetry of both the operators ∂_t and H, so iΓSΓ is a self-adjoint operator (Murphy et al., 2017b), where $i=\sqrt{-\mathrm{1}}$ is the imaginary unit, and $\mathrm{\Gamma }=-\mathbf{\nabla }(-\mathrm{\Delta }{)}^{-\mathrm{1}}\mathbf{\nabla }\cdot$ is the same projection operator arising in the setting of two-component composites. Equation (28) shows that brine advection enhances the thermal diffusivity since ${\mathit{\kappa }}^{*}\ge \mathit{\epsilon }$.

Figure 8 Spectral measures for effective diffusivities. (a) Streamlines for BC flow with fluid velocity field $\mathit{v}=(C\cos y,B\cos x)$ with B=1, $C=\mathrm{1}-p$, and increasing values of p from left to right. BC cell flow with closed streamlines arises for $B=C=\mathrm{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 ε=0. BC shear flow is attained when B=1 and C=0, for which the spectral measure is known to be a δ function at the origin λ=0 (Avellaneda and Majda, 1991), as shown in the rightmost panel of (a) and (c). (b) Corresponding spectral functions ν₂₂(λ) (histogram representations) for the spectral measures ν₂₂ in (c), which display the spectral measure weights $|m_j{|}^{\mathrm{2}}$ vs. eigenvalues λ_j of the matrix M=iΓHΓ for the value of p shown in (a). Since the streamlines do not become elongated in the x direction, the spectral functions for the measure ν₁₁ are qualitatively similar to those of ν₂₂ for BC cell flow, for all p. The spectral functions in (b) are ensemble averaged for B=1 and $C=\mathrm{1}-\mathit{\zeta }$ with $\mathit{\zeta }\sim U\left(\right[\mathrm{0},p\left]\right)$, with p=0.01, 0.5, 0.7, and 1, from left to right. The zoomed-in insets in (c) with $-\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}\le {\mathit{\lambda }}_i\le \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ show the density of the measure near the spectral origin λ=0 increasing as BC cell flow transitions to BC shear flow.

Analytical calculations of the spectral measure ν are extremely difficult except for simple flows like shear flow (Avellaneda and Majda, 1991). However, Padé approximants $[L/M]$ provide rigorous, converging upper and lower bounds (Baker and Graves-Morris, 1996) for the Stieltjes function $f\left(z\right)=\langle |\mathbf{\nabla }{\mathit{\psi }}_k{|}^{\mathrm{2}}\rangle /z=-F/z$ in Eqs. (28) and (29), where $z=\mathrm{1}/{\mathit{\epsilon }}^{\mathrm{2}}$; Padé approximants can be calculated using the moments ν_n of ν. Despite this, the lack of a method to calculate the moments ν_n of ν has impeded progress on obtaining explicit bounds for specific flows using this procedure (Avellaneda and Majda, 1991, 1989) since 1991. We have recently developed a mathematical framework that can be used to compute, in principle, all of the moments ν_n for a spatially or space–time periodic velocity field v, hence Padé 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.

Spectral measure computations for advection diffusion processes

We have extended our numerical methods discussed for two-component media in Sect. 3.1 to compute the spectral measure ν for spatially periodic flows (Murphy et al., 2020a) and developed Fourier methods for computing ν for space-time periodic flows (Murphy et al., 2017b). For simplicity, we focus our discussion to the setting of spatially periodic flows. Computing the spectral measure ν for a given flow involves discretizing the spatially dependent stream matrix H(x), which becomes a banded antisymmetric matrix satisfying ${\mathbf{H}}^T=-\mathbf{H}$. The projection matrix Γ is given by that in Sect. 3.1, and the key self-adjoint operator is given by G=iΓHΓ, which becomes a Hermitian matrix M. In this case, the integral in Eq. (28) and the resolvent in Eq. (30) are given in terms of the eigenvalues λ_j and eigenvectors w_j of the matrix M

$$\begin{array}{}\text{(31)}& \begin{array}{rl}& \mathbf{\nabla }{\mathit{\psi }}_k=\sum _j{\displaystyle \frac{m_j}{\mathit{\epsilon }-i{\mathit{\lambda }}_j}}{\mathit{w}}_j,\langle |\mathbf{\nabla }{\mathit{\psi }}_k{|}^{\mathrm{2}}\rangle =\sum _j〈{\displaystyle \frac{|m_j{|}^{\mathrm{2}}}{{\mathit{\epsilon }}^{\mathrm{2}}+{\mathit{\lambda }}_j^{\mathrm{2}}}}〉,\\ & m_j={\mathit{w}}_j\cdot {\mathit{g}}_j,\end{array}\end{array}$$

which is analogous to Eq. (9).

The computations in Murphy et al. (2017b, 2020a) and those displayed in Fig. 8 show that the spectral origin λ=0 for advection diffusion plays the role of the spectral endpoints $\mathit{\lambda }=\mathrm{0},\mathrm{1}$ for two-component and polycrystalline composite materials. An increase in spectral mass at λ=0 increases the advection-driven enhancement of effective diffusivity above the bare molecular diffusivity ε in the advection-dominated regime, where ε≪1 (or Péclet number ξ≫1). For example, the closed streamlines shown in the leftmost panel of Fig. 8a for BC cell flow, with fluid velocity field $\mathit{v}=(C\cos y,B\cos x)$ and $B=C=\mathrm{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 ${\mathit{\kappa }}^{*}\sim {\mathit{\epsilon }}^{\mathrm{1}/\mathrm{2}}$ (Fannjiang and Papanicolaou, 1994, 1997; Murphy et al., 2020a), vanishing as ε→0. As shown in Fig. 8b and c, and also in Murphy et al. (2017b, 2020a), this is reflected in the spectral measure ν by the lack of adequate mass near λ=0 for the singular integrand $\mathrm{1}/({\mathit{\epsilon }}^{\mathrm{2}}+{\mathit{\lambda }}^{\mathrm{2}})$ to overcome the multiplicative factor of ε for ${\mathit{\kappa }}^{*}=\mathit{\epsilon }(\mathrm{1}+\langle |\mathbf{\nabla }{\mathit{\psi }}_k{|}^{\mathrm{2}}\rangle )$ in Eq. (28).

On the other hand, when B≠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 λ=0 for the singular integrand $\mathrm{1}/({\mathit{\epsilon }}^{\mathrm{2}}+{\mathit{\lambda }}^{\mathrm{2}})$ to overcome the multiplicative factor of ε for ${\mathit{\kappa }}^{*}=\mathit{\epsilon }(\mathrm{1}+\langle |\mathbf{\nabla }{\mathit{\psi }}_k{|}^{\mathrm{2}}\rangle )$ in Eq. (28), leading to a non-zero value of κ^* in the limit ε→0 (Murphy et al., 2017b). This is a key example of how the behavior of the spectral measure ν governs the behavior of the bulk transport coefficient κ^*.

Figure 9 Eigenvalue spacing statistics for sea ice melt ponds (a) and long-range eigenvalue statistics for the brine structures in sea ice (b). (a) Eigenvalue spacing distribution (ESD) P(z) for the melt ponds shown in Fig. 4. (b) Spectral statistics for the brine structures shown in Fig. 3. We see the transition toward universal Wigner–Dyson statistics as melt water phases and brine phases become connected over the scale of the sample.

7 Random matrix theory in sea ice physics

In random matrix theory (RMT) (Guhr et al., 1998; Bohigas and Giannoni, 1984; Deift and Gioev, 2009), long- and short-range correlations of the bulk eigenvalues away from the spectral edge (Canali, 1996; Guhr et al., 1998) for random matrices are measured using various eigenvalue statistics (Guhr et al., 1998; Bohigas and Giannoni, 1984), such as the eigenvalue spacing distribution (ESD), spectral rigidity Δ₃, and number variance Σ². To observe statistical fluctuations of these bulk eigenvalues about the mean density, the eigenvalues must be unfolded (Bohigas and Giannoni, 1984; Guhr et al., 1998; Canali, 1996; Plerou et al., 2002). The localization properties of the eigenvectors are measured in terms of quantities such as the inverse participation ratio (IPR) (Plerou et al., 2002; Evers and Mirlin, 2008).

Figure 10 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 Im(s)=0.001 with $\mathrm{0}\le \mathit{Re}\left(s\right)\le \mathrm{1}$. The electric fields are displayed for values of 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 w_j plotted versus index j. The vertical lines define the δ components of μ, where the eigenvalues satisfy ${\mathit{\lambda }}_j\lesssim {\mathrm{10}}^{-\mathrm{14}}$ and $\mathrm{1}-{\mathit{\lambda }}_j\lesssim {\mathrm{10}}^{-\mathrm{14}}$. The horizontal lines mark the IPR value $I_{\mathrm{GOE}}=\mathrm{3}/N_{\mathrm{1}}$ for the Gaussian orthogonal ensemble (GOE) with matrix size N₁≈ϕN, where N=L^dd.

In Murphy et al. (2017a), we found that 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). The eigenvectors de-localize, and mobility edges appear (Murphy et al., 2017a), similar to the metal/insulator transition in solid state physics. We explored the transition in the 2D and 3D RRN, as well as in sea ice microstructures such as in 2D discretizations of the brine microstructure of sea ice (Golden et al., 1998a, 2007; Golden, 2009), melt ponds on Arctic sea ice (Hohenegger et al., 2012), the sea ice pack itself, and porous human bone (Golden et al., 2011; Kabel et al., 1999; Bonifasi-Lista and Cherkaev, 2009; Cherkaev and Bonifasi-Lista, 2011). We extended these results to two-component composite media with quasiperiodic microgeometry in Morison et al. (2022).

For highly correlated WD spectra exhibited by, for example, real-symmetric matrices of the GOE, the nearest-neighbor ESD P(z) is accurately approximated by $P\left(z\right)\approx (\mathit{\pi }z/\mathrm{2})\exp (-\mathit{\pi }{z}^{\mathrm{2}}/\mathrm{2})$, which illustrates eigenvalue repulsion and vanishes linearly as spacings z→0 (Guhr et al., 1998; Stone et al., 1991; Canali, 1996). In contrast, the ESD for uncorrelated Poisson spectra, $P\left(z\right)=\exp (-z)$, allows for significant level degeneracy (Guhr et al., 1998). In Fig. 9a we display the ESDs for Poisson (blue) and WD (green) spectra, along with the behavior of the ESDs for the matrix M=χ₁Γχ₁, corresponding to the Arctic melt ponds in Fig. 4 with fluid area fraction ϕ. It shows that for sparsely connected systems, the behavior of the ESDs is well described by weakly correlated Poisson-like statistics (Canali, 1996). With increasing connectedness, the ESDs transition toward highly correlated WD statistics with strong level repulsion. This behavior of the ESD reveals a mechanism for the collapse in the spectral gaps of μ. For sparsely connected systems, the weak level repulsion allows for significant level degeneracy and resonances in μ as shown in Murphy et al. (2015) for the 2D percolation model and in Fig. 4a 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 δ components of the measure at the spectral endpoints $\mathit{\lambda }=\mathrm{0},\mathrm{1}$ (as shown in the figures in Sect. 3.1). Our computations of Δ₃ and Σ² are shown in Fig. 9b for the brine microstructure in Fig. 3. Notice that as the system becomes increasingly connected, the long-range statistics transition towards that of the GOE, indicating an increase in the long-range correlations of the eigenvalues.

The eigenvectors w_j of M=χ₁Γχ₁, associated with the N₁×N₁ submatrices of Γ, also exhibit a connectedness-driven transition in their localization properties. The IPR is defined as $I_j={\sum }_i|w_j^i{|}^{\mathrm{4}}$, $i,j=\mathrm{1},\mathrm{\dots },N_{\mathrm{1}}$, where $w_j^i$ is the ith component of w_j. Eigenvectors of matrices in the GOE are known to be highly extended (Deift and Gioev, 2009), with an asymptotic value of the IPR given by $I_{\mathrm{GOE}}=\mathrm{3}/N_{\mathrm{1}}$ (Plerou et al., 2002).

In Murphy et al. (2017a), we found for the 2D and 3D percolation models that as p surpasses p_c and long-range order is established in a RRN, “mobility edges” form in the eigenvector IPR, with a sudden increase in the number of extended eigenvectors. This is analogous to Anderson localization, where mobility edges mark the characteristic energies of the metal–insulator transition (MIT) (Guhr et al., 1998). Remarkably, the mobility edges for RRN are due to very extended eigenstates associated with δ components that form at the spectral endpoints precisely at the percolation threshold p_c (and 1−p_c for 3D) (Murphy and Golden, 2012), which control critical behavior in insulator/conductor and conductor/superconductor systems (Murphy and Golden, 2012; Clerc et al., 1990; Bergman and Stroud, 1992). This and other eigenvector phenomena were observed for two-component composite media with quasiperiodic microgeometry in Morison et al. (2022).

The IPR phenomenon for the eigenvectors of the matrix M=χ₁Γχ₁ is shown for sea ice brine microstructure in Fig. 10c. The electric field E within sea ice brine microstructure exhibits a frequency-dependent Anderson localization transition, as shown in Fig. 10a and b. To generate these figures, the IPR of the electric field, IPR[E](s) in the brine phase χ₁E, was calculated as a function of s via Eq. (9) for Im s=0.001 and $\mathrm{0}\le \mathit{Re}s\le \mathrm{1}$ and 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 Fig. 10a 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.

8 Conclusions

We have gone on a tour of various areas 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, akin to 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, viewing the behavior of these spectral representations through the lens of random matrix theory has revealed fascinating parallels between quantum transport in disordered media, Anderson localization, and classical transport phenomena.