## Abstract

An effective property of a composite material consisting of inclusions within a host matrix depends on the geometry and connectedness of the inclusions. This dependence may be quite strong if the constituents have highly contrasting properties. Here, we consider the inverse problem of using effective property data to obtain information on the geometry of the microstructure. While previous work has been devoted to recovering the volume fractions of the constituents, our focus is on their connectedness—a key feature in critical behaviour and phase transitions. We solve exactly a reduced inverse spectral problem by bounding the volume fraction of the constituents, an inclusion separation parameter and the spectral gap of a self-adjoint operator that depends on the geometry of the composite. We present a new algorithm based on the Möbius transformation structure of the forward bounds whose output is a set of algebraic curves in parameter space bounding regions of admissible parameter values. These results advance the development of techniques for characterizing the microstructure of composite materials. As an example, we obtain inverse bounds on the volume fraction and separation of the brine inclusions in sea ice from measurements of its effective complex permittivity.

## 1. Introduction

We consider the problem of recovering information on the microgeometry of heterogeneous two-phase media from measurements of its effective electromagnetic properties. In the case of highly contrasting phases, as addressed herein, the effective property can depend quite strongly on the connectedness of one of the phases. This may be seen, for example, in the case that the effective property is effective conductivity, and one phase is a good conductor and the other is a good insulator.

Our analysis of this inverse problem is intimately connected with the theory of forward bounds on the effective properties of composite media. Landmarks in the forward theory include the arithmetic and harmonic mean bounds (Wiener 1912), the bounds of Hashin & Shtrikman (1962) and the improved translation bounds (Cherkaev 2000; Milton 2002). The classical derivation of these bounds relies on energy variational principles that are not readily extendible to the interaction of composites with wave fields. If the wavelength of the electromagnetic field is long compared with the scale of the microstructure, then material parameters can assume complex values. Cherkaev & Gibiansky (1994) extended the variationally derived bounds to complex parameters. See also Cherkaev (2000) and Milton (2002).

An alternative approach pioneered by Bergman (1980) and Milton (1980) and developed further by Golden & Papanicolaou (1983) uses analytic continuation for representing and bounding effective properties in the complex case. Analyticity of *m*(*h*), where *m*=*σ**/*σ*_{2} and *h*=*σ*_{1}/*σ*_{2}, is exploited to obtain a Stieltjes integral representation for *F*(*s*)=1−*m*(*h*), where *s*=1/(1−*h*). The integral representation involves the spectral measure *μ* of a self-adjoint operator that depends only on the composite geometry. The support of *μ* lies in the interval [0,1]. The mass of *μ* is *p*_{1}, the volume fraction of medium 1. The analytic continuation method yields a sequence of increasingly tighter bounds that include the arithmetic and harmonic mean bounds and the Hashin–Shtrikman bounds. Information on the geometry of the medium is incorporated into the bounds through the moments of *μ*.

However, these bounds do not incorporate key information about the separation of inclusions, which is important in estimating the effective properties of high-contrast composites. Bruno (1991) advanced the incorporation of separation information within the analytic continuation approach by introducing the class of *matrix particle composites*, and related the separation of the phases to gaps in the support of the spectral measure *μ* at the ends of the interval [0,1].

A matrix particle composite consists of a matrix of conductivity *σ*_{2} containing non-touching inclusions of conductivity *σ*_{1}. A particular case of such a composite is a q-material: spheres or circles of radius *r*_{1} filled with material 1 are surrounded by a corona with outer radius *r*_{2} filled with material 2 and q=*r*_{1}/*r*_{2}<1. As , the corona vanishes, and the particles of material 1 are allowed to touch.

For a q-material in finite rectilinear regions in dimensions 2 and 3, Bruno found explicit formulae for the endpoints of the interval containing the zeros and poles of *m*(*h*), corresponding to the interval of support of the spectral measure *μ*. Bruno applied these results with the analytic continuation method to obtain bounds on the effective conductivity of random matrix particle composites. The bounds yield good estimates on *σ**, even in the limiting high-contrast cases of *h*=0 or , and represent improvements on the classical arithmetic and harmonic mean bounds and the Hashin–Shtrikman bounds. The matrix particle bounds were subsequently extended to media with complex parameters by Sawicz & Golden (1995) and Golden (1998) and applied to estimating the effective complex permittivity of sea ice, extending earlier analysis that made use of the complex elementary and Hashin–Shtrikman bounds.

The inverse problem of recovering information about the microgeometry of a composite from measurements of its effective complex permittivity was first introduced by McPhedran *et al.* (1982) and McPhedran & Milton (1990) in the context of estimating the volume fractions in a two-phase mixture. An analytic approach to bounding microstructural parameters, in particular volume fractions, was developed in Cherkaeva & Tripp (1996), Tripp *et al.* (1998) and Cherkaeva & Golden (1998) and applied to the estimation of brine volume in sea ice in Cherkaeva & Golden (1998) and Gully *et al.* (2007).

The problem is a particular case of *inverse homogenization*: recovery of structural information through identification of the spectral measure *μ* from effective property data. Uniqueness of the spectral measure was established in (Cherkaev 2001) and the theory was further developed in Cherkaev (2003), Cherkaev & Zhang (2003), Cherkaev & Ou (2008), Bonifasi-Lista & Cherkaev (2008) and Zhang & Cherkaev (2009).

Here, we consider a *reduced inverse spectral problem* for matrix particle composites: recovery of bounds on the spectral gap, and hence bounds on q. In particular, we consider q-materials whose high-contrast constituents have volume fractions *p*_{1}=p and *p*_{2}=1−p. Given data on the effective complex permittivity *ε**, the complex matrix particle bounds constructed by Sawicz & Golden (1995) and Golden (1998) are inverted to obtain admissible regions in (p,q)-parameter space, with explicitly computed algebraic curves forming the boundaries. For fixed p, the bounds on q represent the first rigourous inversion for separation or connectedness information on the inclusions in a composite material.

We apply the results to obtain information about the microstructure of sea ice from its effective complex permittivity. One motivation for determining such inverse bounds comes from remote electromagnetic sensing of polar ice packs: brine inclusion separation may be used to monitor the brine percolation threshold (Golden *et al.* 1998*a*, 2007), an *on–off switch* for fluid flow through sea ice that mediates a broad range of biological and climatic processes.

Although we consider sea ice as an example, the inversion algorithm we present may be generalized to other composites and other effective parameters. The simplicity and the generality of the suggested method indicate wide potential applicability and advances in the theory of inverse homogenization for recovering structural parameters from bulk property measurements.

Before embarking on its applications to sea ice, let us describe the general algorithm. The forward bounds of the analytic continuation method all belong to a special class of Möbius transformations. These are complex functions of the form *T*(*z*)=(*Az*+*B*)/(*Cz*+*D*), whose coefficients *A*, *B*, *C*, *D* are usually complex numbers. However, for the Möbius transformations arising in the theory, these coefficients are generally polynomial functions of certain structural parameters—such as p and q—that we are interested in estimating. The problem is to find the ‘good’ parameter values that allow for the forward bounds to capture an observed effective permittivity. This problem is essentially equivalent to a search for those parameter values that allow *ε**=*T*(*z*) to be solved for some real number , known as the spectral parameter. Here, *T*(*z*) describes one of the forward bounding circles. The process of solving *ε**=*T*(*z*) yields an algebraic equation in the structural parameters whose solution bounds the region of ‘good’ parameter values within a larger parameter space. We present the general form of this algebraic equation. In different applications, this algebraic equation will take different forms. If there are two parameters of interest, the solution of this algebraic equation is a one-dimensional curve that lends itself well to visual display.

## 2. Sea ice and its effective complex permittivity

Sea ice is a mixture of pure ice, brine and air. In cold first year ice, brine is typically concentrated in sub-millimetre scale pockets. Although the local complex permittivity *ε*(*x*) varies considerably on the millimetre scale, an electromagnetic wave of sufficiently long wavelength cannot resolve the fine-scale variation between ice, brine and air. This is the case for the Synthetic Aperture Radar used in remote sensing that operates in the C-band—a nominal frequency range 8–4 GHz with a corresponding wavelength range of 3.75–7.5 cm. For such long wavelengths, sea ice may be treated as if it were a homogeneous material having a single effective complex permittivity *ε**. In the long wavelength regime, the wavelength is assumed to be infinite—currently a necessary assumption for the analytic continuation method.

Here, we treat sea ice as a two-component medium instead of a three-component medium. The first phase consists of brine, retained as a pure medium, having complex permittivity *ε*_{1}. The second phase is an ice–air composite having an *effective* complex permittivity *ε*_{2} that is approximated by a mixing formula incorporating the relatively close *ε*_{air} and *ε*_{ice}. Section 6 elaborates on this mixing formula. In the sequel, the second ‘ice phase’ refers to this composite of pure ice and air. Although the three-component case can also be treated with multi-component bounds (Golden 1986; Milton 1987*a*,*b*; Milton & Golden 1990), the mathematics involves several complex variables and has a number of unresolved issues.

The tightest of the forward bounds are obtained in §3*b*(iii) by assuming that the sea ice is both statistically isotropic and is a matrix particle composite. In actual sea ice, the brine inclusions tend to be elongated in the vertical direction. Because of this, we take the dimension *d*=2, so that an assumption of statistical isotropy effectively reduces to an assumption that the geometry is isotropic only within the horizontal plane. This is consistent with the dataset analysed in §7, which comes from measurements of ice slabs that are using vertically incident waves. Except where dimension *d* is explicitly mentioned, all forward and inverse bounds assume *d*=2.

Although this two-dimensional, two-component model works well for sea ice, it may not be viable for other composites such as the Ag–MgF_{2} cermet films studied by Gajdardziska-Josifovska *et al.* (1989), who point out that three-phase cermets are evidently not amenable to accurate modelling by two-phase systems.

## 3. The forward theory

This section summarizes the analytic continuation method (Bergman 1978; Milton 1980; Golden & Papanicolaou 1983; Cherkaev 2001; Zhang & Cherkaev 2009). We also summarize the work of Bruno (1991). We obtain four types of forward bounds on effective complex permittivity: the first- and second-order bounds and , and the first- and second-order matrix particle bounds and . Each of these incorporates different assumptions about the sea-ice geometry: assumes that the brine volume fraction is known; assumes in addition that the distribution of the brine is statistically isotopic. The same holds for and , with the yet additional assumption of a matrix particle model of sea ice.

The analytic continuation method models sea ice as a two-phase random medium in all of **R**^{d}, with an isotropic local complex permittivity *ε*(*x*,*β*), with *ε*(*x*,*β*) a stationary random field in *x*∈**R**^{d} and *β*∈*Ω*, where *Ω* is the set of all realizations of the random medium. The complex permittivity *ε*(*x*,*β*) takes values *ε*_{1} and *ε*_{2}, the permittivities of brine and ice, respectively, and we write *ε*(*x*,*β*)=*ε*_{1}*χ*_{1}(*x*,*β*)+*ε*_{2}*χ*_{2}(*x*,*β*), where *χ*_{1} is the characteristic function of medium 1, which equals one for all realizations *β*∈*Ω* having medium 1 at *x*, and equals zero otherwise, and *χ*_{2}=1−*χ*_{1}.

The constitutive relation can be written as *D*=*εE*, where *E*(*x*,*β*) and *D*(*x*,*β*) are the stationary random electric and displacement fields, satisfying the equations
3.1
We assume that 〈*E*(*x*,*β*)〉=*e*_{j}, where *e*_{j} is a unit vector in the *j*th direction, for some *j*=1,…,*d*, and 〈⋅〉 means ensemble average over *Ω* or spatial average over all of **R**^{d}. The effective complex permittivity tensor ** ε*** is defined as
3.2
Here, we are dealing with isotropic composites, so we consider only one diagonal coefficient of the effective permittivity tensor

*ε**=

*ε**

_{jj}.

Homogeneity of the effective parameter, *ε**(*aε*_{1},*aε*_{2})=*aε**(*ε*_{1},*ε*_{2}), applied with *a*=1/*ε*_{2}, results in *ε**(*ε*_{1}/*ε*_{2},1)=*ε**(*ε*_{1},*ε*_{2})/*ε*_{2}. Hence, by introducing *h*=*ε*_{1}/*ε*_{2}, we can consider the effective complex permittivity formed by constituents with the parameters *h* and 1. We define *m*(*h*) by *m*(*h*)=*ε**(*ε*_{1},*ε*_{2})/*ε*_{2}.

The function *m*(*h*) is analytic off the negative real axis in the *h*-plane and maps the upper half plane to the upper half plane: this characterizes *m*(*h*) as a Herglotz function (Baker & Graves-Morris 1996, p. 262). A Herglotz function has the integral representation
with *ν*(*u*) bounded and non-decreasing on **R**. If the first moment of *ν* is finite, we may rewrite this as

This representation allows us to write an analytic integral representation for *ε**. For this purpose, it is more convenient to introduce *s*=1/(1−*h*) and consider *F*(*s*)=1−*m*(*h*), which is analytic off [0,1] in the *s*-plane. Then,
3.3
where *μ* is a positive measure on [0,1]. This formula is essentially the spectral representation of the resolvent *E*=(*s*+*Γχ*_{1})^{−1}*e*_{j}, obtained from (3.1) and (3.2), where *Γ*=∇(−Δ)^{−1}∇⋅, Δ denotes the Laplacian, and (−Δ)^{−1} is convolution with the free-space Green function for −Δ. In the Hilbert space *L*^{2}(*Ω*,*P*) with weight *χ*_{1} in the inner product, *Γχ*_{1} is a bounded self-adjoint operator with norm less than or equal to one. In (3.3), *μ* is a spectral measure of *Γχ*_{1}. One of the most important features of (3.3) is that it separates the parameter information in *s*=*ε*_{2}/(*ε*_{2}−*ε*_{1}) from information about the geometry of the mixture, which is all contained in *μ*.

Statistical assumptions about the geometry are incorporated into *μ* through its moments *μ*_{n}. A comparison of the perturbation expansion of (3.3) around a homogeneous medium, or equivalently *ε*_{1}=*ε*_{2}, namely
3.4
with a similar expansion of the resolvent representation for *F*(*s*) yields . The zeroth moment *μ*_{0}=*p*_{1}, where *p*_{1} is the volume fraction of the first material in the composite, and for a statistically isotropic composite material, the first moment of *μ* is calculated as *μ*_{1}=*p*_{1}*p*_{2}*d*^{−1}.

Expansion (3.4) converges only in the disc |*h*−1|<1, while the integral representation (3.3) provides the analytic continuation of (3.4) to the full domain of analyticity. In this way, information obtained about a nearly homogeneous system can be used to analyse the system near percolation as or .

### (a) First- and second-order bounds

Bounds on *ε**, or *F*(*s*), are obtained by fixing *s* in (3.3), varying over admissible measures *μ*, and finding the corresponding range of values of *F*(*s*) in the complex plane. Two types of bounds on *ε** are readily obtained. The first-order bounds assume only that the relative volume fractions *p*_{1} and *p*_{2}=1−*p*_{1} are known, so that only *μ*_{0}=*p*_{1} need to be satisfied. The second-order bounds assume in addition that the material is statistically isotropic.

#### (i) The forward region

The set of admissible measures satisfying *μ*_{0}=*p*_{1} forms a compact, convex set. Since (3.3) is a linear functional of *μ*, the extreme values of *F* are attained by extremes of the admissible measures: these are the Dirac point measures of the form *p*_{1}*δ*_{z}. The values of *F*(*s*) are, therefore, bounded by the circle
3.5
A second circle bounding may be obtained by introducing the function
Bergman (1982) established that this is a Herglotz function like *F*(*s*), analytic off [0,1], with a representation like (3.3) whose representing measure has mass *p*_{2},
3.6
Then, in the *E*-plane, the other circular boundary of is parametrized by
3.7
After these circles are transferred to the common *ε**-plane, and noting their intersection, the relevant bounding arcs given parenthetically in (3.5) and (3.7) may be determined.

#### (ii) The forward region

If the material is further assumed to be statistically isotropic, i.e. *ε**_{ik}=*ε***δ*_{ik}, then *μ*_{1}=*p*_{1}*p*_{2}*d*^{−1} must be satisfied as well, and we obtain the second-order bounds . Now instead of directly varying over all those admissible measures whose moments satisfy both *μ*_{0}=*p*_{1} and *μ*_{1}=*p*_{1}*p*_{2}*d*^{−1}, it is convenient to use the following transformation, as pointed out by Bergman (1982):
3.8
The function *F*_{1}(*s*) is also an Herglotz function, and has the representation
3.9
The constraints on the moments *μ*_{0} and *μ*_{1} are then transformed to a restriction only on the mass of *μ*^{1}, which is . Applying the same procedure that was used for yields the region , whose boundaries are again circular arcs. One arc, in the *F*-plane, is
3.10
In the *E*-plane, the other arc is
3.11

The forward bounds discussed up to this point are summarized in table 1: on the left-hand side are the bounds on the values *F*(*s*) and *E*(*s*), which arise from the integral representations of the form (3.3) and (3.6) by letting the representing measure vary over the extremes. In order to be useful as bounds on the complex permittivity, these should be transferred to the *ε**-plane using the relations
3.12
The functions *T*_{i,j}(*z*;p), *i*,*j*=1,2, defined on the right-hand side of table 1 are obtained by transferring such bounds accordingly. Some notational simplification is achieved by introducing the new variable *θ*=2*s*−1. Note that each *T*_{i,j}(*z*;p) is a Möbius transformation that maps the extended real line onto a circle in the *ε**-plane.

#### (iii) Isomorphic groups

A convenient way of transferring bounds between planes exploits the isomorphism between the Möbius transformation group and the projective general linear group PGL(2,**C**) (Jones & Singerman 1987). For example, the relationship between *ε** and *F* given by (3.12) has the matrix representations, with *s* fixed,
Then, with p=*p*_{1}, multiplying a matrix (in the equivalence class) representing *C*_{1}(*z*) on the left by yields a matrix representation for *T*_{1,1}(*z*;p),
The other entries on the right-hand side of table 1 may be obtained similarly. Let us note that and .

#### (iv) The classical bounds

In the *ε**-plane, the two arcs describing meet at the vertices
and
When *ε*_{1} and *ε*_{2} are real, reduces to the interval defined by the classical harmonic mean and arithmetic mean bounds: (*p*_{1}/*ε*_{1}+*p*_{2}/*ε*_{2})^{−1}≤*ε**≤*p*_{1}*ε*_{1}+*p*_{2}*ε*_{2}. Similarly, when *ε*_{1} and *ε*_{2} are real with *ε*_{1}≥*ε*_{2}, collapses to the interval
the very bounds derived by Hashin & Shtrikman (1962).

### (b) First- and second-order matrix particle bounds

Still tighter bounds on *ε**, and may be obtained if the material has a matrix particle structure with separated inclusions. In this case, the support of the spectral measure *μ* in (3.3) is contained in the subinterval [*s*_{m},*s*_{M}] and
3.13
This important observation is due to Bruno (1991), whose work we summarize next. Although Bruno considers the effective conductivity of strongly heterogeneous composites, the mathematics of effective complex permittivity is similar.

#### (i) Summary of Bruno's work

The bounds *s*_{m} and *s*_{M} are obtained by mapping the bounds on the singularities and zeros of the corresponding function *m*(*h*) from the complex *h*-plane to the *s*-plane. The *h*-plane bounds result from analyticity of *m*(*h*), corresponding to a matrix particle composite, in neighbourhoods of *h*=0 and . The method is based on an expansion of the electric potential *ϕ* into convergent series: in *h* around *h*=0; in *w*=1/*h* around .

To sketch the argument, we consider a domain large in comparison with the size of inhomogeneity, , filled with a composite mixture of two materials with properties *ε*_{1}=*h* and *ε*_{2}=1. We assume that an electric potential *ϕ* is constant on the upper ∂_{1} and lower ∂_{0} boundaries and periodic on vertical boundaries *S* of the domain , so that *ϕ* satisfies the equation
3.14
Using the solution of (3.14), the effective conductivity *m*(*h*) may be calculated as
3.15
where is the part of occupied by the material with conductivity *h*, and . It is known that (3.14) has a unique solution for any value of *h* outside the negative real axis (Bergman 1978; Golden & Papanicolaou 1983) and that *m*(*h*) is an analytic function of .

Using a combination of the trace and extension theorems (Nečas 1967) and the Poincaré inequality for functions in and , it can be shown that for the matrix particle composites, the two integrals on the right-hand side of (3.15) are bounded: each function can be extended to a function , which coincides with *u*′ in and
3.16
and for each function , there is a function that differs by a constant from *u* on the boundary of each inclusion, and
3.17
The fields in matrix particle composites with coronas around grain inclusions satisfy these bounds with positive constants *A* and *B* depending on grain shape and separation.

To show analyticity of *m*(*h*) around the origin, the solution of (3.14) is represented as
3.18
with the coefficients *ϕ*_{k} satisfying an infinite system of equations obtained from (3.14) upon substitution of the series (3.18). Separating the problem into two problems over subdomains and allows one to solve them iteratively, since the problems are coupled through the condition on the boundary of inclusions. Using (3.16), one can show that (3.18) as well as the series of the gradients ∇*ϕ*_{k} is absolutely convergent for |*h*|≤1/*A*. Therefore, *ϕ*(*x*,*h*) and hence *m*(*h*) are analytic for |*h*|≤1/*A*. Similarly, to show analyticity of *m*(*h*) at infinity, *ϕ*(*x*,*h*)=*ψ*(*x*,*w*=*h*^{−1}) is expanded as
3.19
around the origin in the *w*-plane. Using (3.17), the convergence of (3.19) can be shown for |*w*|≤1/*B*.

Together, (3.18) and (3.19) provide the analytic continuation of *ϕ*(*x*,*h*) in the regions |*h*|≤1/*A* and |*h*|≥*B*. Accordingly, *m*(*h*) is analytic in these regions, and in particular, for *h* on the negative real axis, −1/*A*≤*h* and *h*≤−*B*. Mapping the interval [−*B*,−1/*A*] to the *s*-plane gives the subinterval *s*_{m}≤*z*≤*s*_{M} bounding the support of the measure *μ* displayed in (3.13), with *s*_{m} and *s*_{M} depending on the microgeometry of the composite.

The further the separation of the inclusions, the smaller the support interval [*s*_{m},*s*_{M}], and the tighter the bounds. Explicit calculations for *s*_{m} and *s*_{M} are available for the following matrix particle model of sea ice: within the horizontal plane, the brine is assumed to be contained in separated, circular discs of radii *r*_{brine}, holding random positions in such a way that each disc of brine is surrounded by a corona of ice, with outer radius *r*_{ice}. As introduced in §1, this is a q-material with q=*r*_{brine}/*r*_{ice}. For such a geometry, Bruno has calculated for *d*=2 that
3.20
Smaller q values indicate well-separated brine—and colder temperatures as illustrated in §7—and q=1 corresponds to no restriction on the separation, with *s*_{m}=0 and *s*_{M}=1, so that and coincide with and .

Returning to (3.13), a convenient way of incorporating the support restriction is the substitution
3.21
which maps [*s*_{m},*s*_{M}] in the *s*-plane to [0,1] in the *t*-plane. Consequently,
3.22
is analytic off [0,1] in the complex *t*-plane, and there is a positive Borel measure *ν* on [0,1] such that . Using *μ*_{0}=*p*_{1} and *μ*_{1}=*p*_{1}*p*_{2}*d*^{−1}, it can be shown that *ν* has moments
3.23a
and
3.23b
with the formula for *ν*_{1} holding under the assumption of statistical isotropy.

#### (ii) The forward region

The bounds are obtained by assuming (3.23a) is satisfied. Applying the same extremal procedure that gave (3.5) for shows that the values of *H*(*t*) lie inside the circle
Note that *H*(*t*) assumes values in the *F*-plane: the *H*- and *F*-planes coincide. This translates via (3.21) and (3.22) into the circle , *z*∈**R**, which is readily seen to coincide with the image of *C*_{1}(*z*). So in this case, the matrix particle assumption provides no improvement over the first arc of .

Next, we consider the analogue *G*(*t*) of *E*(*s*) defined by
3.24
Then, *G*(*t*) has an integral representation , where the mass of *ρ* is *ρ*_{0}=1−*λ*^{−1}*p*_{1}. We then obtain a circle in the *G*-plane analogous to in (3.7), namely,
3.25
In the *F*-plane, this becomes, via (3.22) and (3.24),
3.26
which improves upon (3.7). Here, (3.26) is a correction of Golden (1998, eqn 3.29).

#### (iii) The forward region

Finally, we consider the case where the material is further assumed to be statistically isotropic. Let
3.27
Then, *H*_{1}(*t*) is a Herglotz function that is analytic off [0,1] and has an integral representation in terms of a measure *ν*^{1}, which can be shown to have mass
3.28
The allowed values of *H*_{1}(*t*) are contained inside the circle , which in the *H*-plane becomes the circle
3.29

The other circle is obtained by applying a similar transformation to *G*(*t*),
3.30
This function is also Herglotz, analytic off [0,1], with representing measure *ρ*^{1}. The mass of *ρ*^{1} is
3.31
The allowed values of *G*_{1}(*t*) are contained inside the circle , which in the *G*-plane becomes
3.32

Table 2 summarizes the first- and second-order matrix particle bounds: the left-hand side gives bounds on the values of *G*(*t*) and *H*(*t*), which are then transferred to the *ε**-plane on the right-hand side. For example, the matrix representation of *T*_{2,2}(*z*;p,q) may be obtained by multiplying the matrix representing on the left by , i.e.
and then substituting *t*=(*θ*+q^{2})/(2q^{2}), *ν*_{0}=p/q^{2} and *ν*_{1}=p(q^{2}−p)/(2q^{4}) from (3.20)–(3.23). Typical forward bounds on *ε** are illustrated in figure 1. While figure 2 illustrates inverse bounds, figure 3 also illustrates forward bounds.

## 4. The inverse problem

Having computed the forward bounds, we now formulate the general inverse problem: given an observed effective property, determine the range of parameter values consistent with the observation. This is illustrated in figure 1 using sea ice modelled as a two-component random medium, effective complex permittivity *ε** as the observed effective property and brine volume fraction p as the parameter of interest. As the brine volume increases from p≈0.007 to p≈0.040, the forward region changes size and sweeps over the fixed observed *ε**, while the smaller region covers *ε** only for p between p≈0.012 and p≈0.024. These two intervals, respectively, are the first- and second-order inverse bounds on p, given the observed *ε**. It turns out that the end points of these intervals may be determined as roots of the polynomials *F*_{i,j}(p) defined in §5 by inverting *T*_{i,j}(*z*;p), *i*,*j*=1,2.

Simultaneous inversion for both p and q is more complicated, but the same method applies: given an observed *ε**, inversion of *T*_{i,j}(*z*;p,q) determines the algebraic curve *G*_{i,j}(p,q)=0 that bounds a region in (p,q)-parameter space.

## 5. The inverse algorithm

Theorem 5.1 provides the theoretical basis of the inverse algorithm. It asserts that the boundary of the inverse region may be computed by solving an algebraic equation involving structural parameters of the measure *μ*. This extends and provides an alternative to the approach of Cherkaeva & Golden (1998), in which bounds for p were found by solving coupled nonlinear equations.

Theorem 5.1 is based on the observation that the bounds on *ε** on the right-hand sides of tables 1 and 2 are all Möbius transformations in *z*; the relevant arcs shown in figure 1 are the images of certain subintervals of [0,1] under such transformations. The full circle is the image of the extended real line . For generality, we now designate structural parameters of the measure *μ* by *π*_{1},…,*π*_{n}, instead of just p and q. If parameter values are chosen so that *ε** (or another observed effective property) lies on such a circle, the spectral parameter is defined to be the such that . We explicitly do not regard as a ‘structural parameter of the measure *μ*’.

We use the following notation: **C**[*x*_{1},…,*x*_{n}] denotes the ring of multi-variate polynomials in *x*_{1},…,*x*_{n}, with complex coefficients; is the complex conjugate of *z*∈**C**; and both *T*_{π}(*z*) and *T*(*z*;*π*) are used to indicate the parametric dependence of a Möbius transformation on *π*=(*π*_{1},…,*π*_{n}). We use the following facts: Möbius transformations map circles to circles in **C**, where lines are regarded as circles of infinite radii and the compositional inverse of *T*(*z*)=(*Az*+*B*)/(*Cz*+*D*) is *T*^{−1}(*z*)=(*Dz*−*B*)/(−*Cz*+*A*) (Jones & Singerman 1987).

### Theorem 5.1

*Let T(z;π) be an n-parameter family of Möbius transformations
*
*where A(π),B(π),C(π),D(π)∈**C**[π*_{1}*,…,π*_{n}*] and the parameters π*_{1}*,…,π*_{n} *are real. Then, a fixed ζ∈**C**lies on* *a circle with a single point removed, if and only if* *is a real root of the multi-variate polynomial F*_{ζ}*(π) defined by
*
5.1
*and* *is not also a root of the multi-variate polynomial
*
5.2
*For such* *, where the spectral parameter* *may be computed as
*
5.3

### Proof.

Let *ζ*∈**C** be fixed and assume , with . Then,
5.4
for some , . Since *S*_{ζ}(*π*)=0 can only occur if , it follows that *S*_{ζ}(*π*)≠0. Solving for by inverting (5.4) gives
5.5
By considering separately the real and imaginary parts of (5.5) and using the fact that , it follows that is a solution of
5.6
In view of , this is equivalent to being a root of the polynomial (5.1).

In the other direction, let *ζ*∈**C** be fixed and suppose and . Then, (5.6) holds with . Let be defined by (5.5). Then, . Since (5.6) holds with , it follows that is real, so (5.5) defines as the spectral parameter. The remaining part of (5.3) holds because for any *z*,*w*∈**C**/{0}, we have if and only if the points {0,*z*,*w*} are collinear, in which case, *z*/*w*=ℜ(*z*)/ℜ(*w*)=ℑ(*z*)/ℑ(*w*). ■

Using this unified approach, we first apply theorem 5.1 to obtain formulae for first- and second-order inverse bounds on the volume fraction p of brine in sea ice, and in passing rederive results of Cherkaeva & Golden (1998). Next, we apply theorem 5.1 to obtain first- and second-order inverse bounds in (p,q)-parameter space. Since the role of *ζ* will always be assumed by *ε**, to simplify notation, we will write *F*(p) for *F*_{ε*}(p) and *G*(p,q) for *F*_{ε*}(p,q).

### (a) Inverting

The forward region is derived assuming only knowledge of the volume fraction p=*p*_{1}. We invert to obtain bounds on p by considering *T*_{1,1}(*z*;p) and *T*_{1,2}(*z*;p) as families of maps parametrized by p.

#### (i) Inverting the first arc of

The polynomial coefficients of *T*_{1,1}(*z*;p) are
An application of theorem 5.1 gives . Solving *F*_{1,1}(p)=0 gives a lower bound on the volume fraction,
5.7
Note that the subscript ‘1’ on refers to its role as a first-order bound, while *p*_{1} denotes a volume fraction.

#### (ii) Inverting the second arc of

The polynomial coefficients of *T*_{1,2}(*z*;p) are
Theorem 5.1 gives , and solving *F*_{1,2}(p)=0 gives an upper bound on the volume fraction,
5.8
The fact that (5.7) gives a lower bound while (5.8) gives an upper bound was determined numerically, i.e. it is the particular values of *ε*_{1}, *ε*_{2}, *ε** that imply .

### (b) Inverting

The forward region is the region in the *ε**-plane bounded by the two circles described by *T*_{2,1}(*z*;p) and *T*_{2,2}(*z*;p). Following the same procedure for as , we find that each circle determines a polynomial that is quadratic in p, yielding the second-order bounds and with .

#### (i) Inverting the first arc of

Inverting *T*_{2,1}(*z*;p) gives *F*_{2,1}(p)=*a*_{1}p^{2}+*b*_{1}p+*c*_{1}, with
Solving *F*_{2,1}(p)=0 gives the second-order upper bound
5.9

#### (ii) Inverting the second arc of

Inverting *T*_{2,2}(*z*,p) gives *F*_{2,2}(p)=*a*_{2}p^{2}+*b*_{2}p+*c*_{2}, with
Solving *F*_{2,2}(p)=0 determines the second-order lower bound
5.10
The fact that (5.9) provides an upper bound while (5.10) provides a lower bound was determined numerically, as was the determination of the sign of the radical in these equations.

### (c) Inverting

The forward region in the *ε**-plane is determined by the two intersecting circles: *T*_{1,1}(*z*;p,q) and *T*_{1,2}(*z*;p,q), *z*∈**R**.

#### (i) Inverting the first arc of

Inversion of *T*_{1,1}(*z*;p,q) reduces to (5.7) because the circular image of coincides with the image of *C*_{1}(*z*), as noted in §3*b*(ii).

#### (ii) Inverting the second arc of

Inverting *T*_{1,2}(*z*;p,q) gives the polynomial
which is cubic in q^{2}. The algebraic curve *G*_{1,2}(p,q)=0 is readily expressed as

### (d) Inverting

The forward region in the *ε**-plane is described by the intersection of the two circles: *T*_{2,1}(*z*;p,q) and *T*_{2,2}(*z*;p,q), *z*∈**R**.

#### (i) Inverting the first arc of

Inverting *T*_{2,1}(*z*;p,q) gives *G*_{2,1}(p,q)=*a*_{1}(q)p^{2}+*b*_{1}(q)p+*c*_{1}(q) with
Note that *G*_{2,1}(p,q)|_{q=1}=*F*_{2,1}(p). The algebraic curve *G*_{2,1}(p,q)=0 is given by
5.11
The curve *G*_{2,1}(p,q)=0 is illustrated in figure 2; the top and bottom portions of the curve correspond to the positive and negative roots in (5.11), respectively.

Alternatively, and this will be used in §7, 5.12

#### (ii) Inverting the second arc of

Inverting *T*_{2,2}(*z*;p,q) gives *G*_{2,2}(p,q)=*a*_{2}(q)p^{2}+*b*_{2}(q)p+*c*_{2}(q), with
The algebraic curve *G*_{2,2}(p,q)=0 is given by
5.13
Both roots appear in the full display of the curve on the right-hand side of figure 2.

### (e) Matching forward and inverse regions

The curves *G*_{i,j}(p,q)=0 are boundary curves; it remains to determine which side of the boundaries the admissible parameter values lie. When q=1, reduces to so that the inverse region determined by must contain the line segment [(p_{1,l},1),(p_{1,u},1)]. Similarly, the inverse region determined by must contain the shorter segment [(p_{2,l},1),(p_{2,u},1)]. A complete matching of the forward regions with the corresponding inverse regions may be carried out by computing forward regions at selected (p,q) pairs. Figure 3 illustrates this process.

### (f) If *F*_{ζ}(*π*)=0 has no real roots

We return to a discussion of the inversion algorithm in general, not just its application to sea ice. We follow the same notation of theorem 5.1 and again consider *π*_{1},…,*π*_{n} as variable (real) parameters.

It may happen that *F*_{ζ}(*π*)=*F*_{ζ}(*π*_{1},…,*π*_{n})=0 has no real roots. The following theorem addresses this issue. Here is its practical application: suppose for an observed effective property *ζ*∈**C**, we have *F*_{ζ}(*π*)≠0 for all relevant parameter values *π*. According to theorem 5.1, this means that the circles *T*_{π}(**R**) do not come into contact with *ζ* as *π* varies. Theorem 5.2 gives criteria to decide if this is because *ζ* always lies inside, or always lies outside, each of the circles *T*_{π}(**R**).

The new notation *Λ*_{n}, denoting a subset of **R**^{n}, accounts for the fact that relevant parameter values may not comprise all of **R**^{n}. For example, for our sea-ice problem, the relevant values of (p,q) lie in . We also introduce the notation denoting a connected subset of *Λ*_{n}. In most applications, we would expect . Jones & Singerman (1987, p. 28) discuss the ‘circle inversion’ mentioned after (5.16).

### Theorem 5.2

*Suppose* *is a connected subset of parameter set Λ*_{n}*⊆**R*^{n} *and
*
5.14
*If F*_{ζ}*(π)=0 has no real roots, then ζ lies outside each circle in the family* *provided the following inequality holds for at least one*
5.15
*If F*_{ζ}*(π)=0 has no real roots, then ζ lies inside each circle in* *, provided that the reverse inequality in (5.15) holds for at least one* *.*

### Proof.

First, the hypothesis (5.14) implies that each *T*_{π}(**R**) is a circle of finite radius. Otherwise, would be a straight line in **C** passing through in the extended complex plane, implying that for some . Upon computing , we find *z*=−*D*(*π*)/*C*(*π*)∉**R** by (5.14), and also by (5.14). Therefore, *T*_{π}(**R**) has finite radius.

Next, let *ξ*(*π*) and *ρ*(*π*) denote the centre and radius of *T*_{π}(**R**), respectively,
5.16
These may be derived by noting that inversion in the circle *T*_{π}(**R**) exchanges *ξ* and , and is conjugate via *T*_{π} to a reflection across **R**. Since , it follows that , yielding the indicated formula; and . Formulae (5.16) show that *ξ*(*π*) and *ρ*(*π*) are continuous. Define by *J*(*π*)=|*ζ*−*ξ*(*π*)|−*ρ*(*π*). Note that . (Otherwise, *ζ*∈*T*_{π}(**R**) for some , and then *F*_{ζ}(*π*)=0 would have a real root in by theorem 5.1.) Since *J* is continuous, is connected. Hence, either *J*(*π*)>0 for all or for all . If (5.15) holds for a particular , then *J*(*π**)>0, hence *J*(*π*)>0 for all , so that *ζ* lies outside each circle *T*_{π}(**R**) for all . Similar reasoning handles the reverse inequality in (5.15). ■

Theorem 5.2 has the following corollary in its specialization to sea ice: if for a given *ε**≠0, *F*_{i,j}(p) fails to have a real root, then *ε** lies *outside* the circles *T*_{i,j}(**R**;p) for all p∈[0,1]; if for a given *ε**≠0, *G*_{i,j}(p,q) fails to have a real root, then *ε** lies *outside* the circle *T*_{i,j}(**R**;p,q) for all (p,q)∈[0,1]×[0,1]. This is because for each entry on the right-hand side of table 1, there exists a p∈[0,1], such that *A*(p)*D*(p)−*B*(p)*C*(p)=0; and for each entry on the right-hand side of table 2, there exists a (p,q)∈[0,1]×[0,1], such that *A*(p,q)*D*(p,q)−*B*(p,q)*C*(p,q)=0.

## 6. Inverse bounds for volume fraction

In this section, we consider the problem of inverting for the single parameter p before considering both p and q in §7. We apply the inversion method developed in §5 to obtain upper and lower bounds on p using the effective complex permittivity data for sea ice given in §8. We then compare these bounds with an empirical formula from Frankenstein & Garner (1967) that determines brine volume from temperature and salinity.

The inversion formulae for , , , depend on *ε**, and the complex permittivities of the two phases, *ε*_{1} and *ε*_{2}. The latter two may be calculated from the measured temperature, sample salinity, sample density and the electromagnetic frequency of the experiments, using the formulae described next.

Concerning *ε*_{1}, the complex permittivity of the brine, we use the calculations of Stogryn & Desargant (1985) that are based on a Debye-type relaxation equation,
6.1
in which *f* denotes the frequency in gigahertz, and and *σ* are expressed as functions of temperature by the following equations fit to experimental data:
Here, *ε*_{s} and are the limiting static and high-frequency values of the real part of *ε*_{1}, *τ* is the relaxation time in nanoseconds, *ε*_{0} is the permittivity of free space, 8.85419×10^{−12} F m^{−1}, and *σ* is the ionic conductivity of the dissolved salts in siemens per metre. It is assumed that *σ* is independent of frequency.

Concerning *ε*_{2}, it was found by Cherkaeva & Golden (1998) that the theoretical forward bounds fit the data more closely by accounting for air in the sea ice. In particular, the complex permittivity of the ice *ε*_{2} was calculated as a permittivity of a composite with a small volume fraction of air using the Maxwell–Garnett formula. We use the same approach here,
6.2
Here, *ε*_{ice} = (3.1884 + 0.00091 *T*) + 0.00005i (Mätzler & Wegmuller 1987, 1988), *ε*_{air}=1, and the volume fraction of air, *p*_{air}, is calculated via the equations given by Cox & Weeks (1983),
Here, *ρ* is the density of the sea-ice sample in grams per cubic centimetre; *T* is its temperature in degrees celsius; *S* is its salinity in parts per thousand; *ρ*_{ice} is the density of pure ice in grams per cubic centimetre, which is given by *ρ*_{ice}=0.917−1.403×10^{−4}*T*; and the coefficients of *F*_{j}(*T*)=*α*_{0}+*α*_{1}*T*+*α*_{2}*T*^{2}+*α*_{3}*T*^{3}, *j*=1,2, are given in table 3. In (6.2), we take *d*=3, as the air inclusions in actual sea ice are uniformly and isotropically distributed throughout the ice in three dimensions, as opposed to the brine inclusions.

Table 4 and figure 4 compare the result of inverting for brine volume fraction from effective complex permittivity with the result of computing brine volume fraction using the equation of Frankenstein & Garner (1967),
6.3
Here, *T* is the temperature in degrees, celsius and *S* is the salinity in parts per thousand.

## 7. Inverse bounds for inclusion separation

Figure 2 shows typical inverse bounds on p and q. For a given brine volume fraction , we may determine an interval of admissible q values from the second-order matrix particle bounds: it is the interval , where is the value of q where a horizontal line at level p will intersect the inverse boundary curve from the first arc of . Thus, may be computed by setting (5.12) equal to zero and solving for q,
7.1
If p is not in the indicated interval, then we set . Here, *d*_{1}(p), *e*_{1}(p) and *f*_{1}(p) are given by (5.12). The need for the minus sign on the square root was established numerically.

Figure 5 shows calculated by (7.1), with p_{c} calculated by (6.3), using data given in table 5. It indicates a coalescence towards percolation as the temperature rises. Since 1 is always the upper bound, we cannot make inferences about the ice being bounded away from percolation, at colder temperatures. Nevertheless, the lower bound is still informative; it bounds q towards percolation.

## 8. Data

Tables 5 and 6 record data used herein; see also the electronic supplementary material.

## Acknowledgements

We gratefully acknowledge support from the Division of Mathematical Sciences and the Office of Polar Programs at the US National Science Foundation through grants DMS-0940249 and ARC-0934721. The authors are grateful to the anonymous referees for useful comments.

- Received September 1, 2011.
- Accepted October 17, 2011.

- This journal is © 2011 The Royal Society