## Abstract

Minimization variational principles for linear elastodynamic, acoustic or electromagnetic time-harmonic waves in dissipative media were obtained by Milton *et al.* (Milton *et al.* 2009 *Proc. R. Soc. A* **465**, 367–396 (doi:10.1098/rspa.2008.0195)), generalizing the quasistatic variational principles of Cherkaev and Gibiansky (Cherkaev & Gibiansky 1994 *J. Math. Phys.* **35**, 127–145 (doi:10.1063/1.530782)). Here, a further generalization is made to allow for a much wider variety of boundary conditions, and in particular Dirichlet and Neumann boundary conditions. In addition minimization or maximization principles of the Hashin–Shtrikman type, incorporating ‘polarization fields’, are developed. The corresponding principles for static problems have found substantial use in bounding the effective static properties of composite materials. The new dynamical principles offer the prospect of developing bounds on the effective dynamic properties of such materials.

## 1. Introduction

Stationary variational principles for wave equations have been the subject of much attention (e.g. Oden & Reddy 1983; Lazzari & Nibbi 2000; Altay & Dökmeci 2004; and references therein). Among these, the variational principles of Gurtin (1964), Tao (1966) and Willis (1981, 1984) minimize variational principles in the Laplace domain but not in the frequency domain. In the frequency domain in a dissipative medium, they correspond to saddle-point variational principles (Borcea 1999; Milton *et al.* 2009).

Saddle-point variational principles for time-harmonic waves in dissipative media in the quasistatic limit, where the wavelengths and attenuation lengths are large compared with the size of the body, were derived by Cherkaev & Gibiansky (1994). They realized that by making a partial Legendre transform, one could convert them to minimization principles. Equivalently, one could rewrite the complex equations in a form involving the real and imaginary parts of the fields with a real positive-definite matrix entering the constitutive law. From this reformulation, and from the differential constraints on the fields, one immediately obtains the minimization principles. Their variational principles were employed to give a simple derivation (Milton 1990) of existing bounds on the effective complex electrical permittivity of lossy multi-phase composites. (These bounds, which generalized the two-phase bounds of Milton (1981) and Bergman (1982), were first conjectured by Golden & Papanicolaou (1985) and Golden (1986) and subsequently proved by Bergman (1986), Milton (1987) and Milton & Golden (1990).) Their variational principle also provided entirely new bounds on the complex bulk and shear moduli of viscoelastic two-phase composites (Gibiansky & Lakes 1993, 1997; Gibiansky & Milton 1993; Gibiansky *et al.* 1999), which, in turn, led to bounds on the complex thermal expansion coefficient of viscoelastic two-phase composites (Berryman 2009).

Milton *et al.* (2009) realized that the technique of Cherkaev and Gibiansky could be extended to obtain minimization variational principles for time-harmonic waves in a body composed of dissipative material. In particular, they obtained minimization variational principles for acoustics, elastodynamics and electromagnetism. (In passing, we mention that the final paragraph in that paper should be disregarded as the inequality (6.8) does not provide any constraints on the moduli inside the body as the difference between the right- and left-hand sides can be expressed as the integral of a square.) An unappealing feature of these variational principles is that they require boundary conditions on the fields that are difficult to physically impose, such as fixing the real part of the displacement together with the real part of the traction (for elastodynamics) or the real part of the tangential component of *E* together with the imaginary part of *H* (for electromagnetism).

Here, by including what amounts to additional surface terms, we derive more general variational principles that allow for all sorts of boundary conditions. These include the standard boundary conditions where both real and imaginary parts of the displacement, traction, tangential component of *E* or tangential component of *H* are specified. Like the variational principles of Milton *et al.* (2009), these variational principles should be useful for tomography where one seeks information about the moduli inside a body from a series of measurements of the fields at the surface of the body from various different boundary conditions.

We remark that the formulation in terms of minimization principles enables one to use the conjugate gradient algorithm for the numerical solution using finite elements. This has been explored by Richins & Dobson (in preparation) for the Helmholtz equation, who also obtained error bounds on the difference between the numerical solution and the exact solution as a function of the grid size.

Another completely different type of variational principle was derived by Hashin & Shtrikman (1962*a*,*b*,*c*, 1963) to bound the effective moduli of composites. This type of principle was extended to wave equations as stationary variational principles by Ben-Amoz (1966) and Willis (1981, 1984). Here, we derive Hashin–Shtrikman type minimization principles for time-harmonic waves in dissipative elastodynamic media. The extension of these principles to acoustics and electrodynamics exactly parallels the treatment given for elastodynamics, and is therefore omitted.

## 2. Preliminaries for elastodynamics

The propagation of stress waves in any body is governed by the equation of motion
2.1
where *σ* denotes the stress tensor, *p* is the momentum density (a vector) and *f* is the body force per unit volume. The body is characterized by its constitutive relations, which will be taken to have the linear form
2.2
where the strain tensor *e* and velocity vector *v* are related to displacement *u* by
2.3
The asterisk denotes convolution with respect to time and thus the body is taken to be viscoelastic. Although it is usually the case that the mass density *ρ* is a real scalar, the formulation is not complicated if *ρ* is taken to be a tensor-valued convolution operator, and therefore this is assumed.

The case of time-harmonic disturbances will be assumed, in which case all fields become functions of the spatial variable *x* only, multiplied by the factor *e*^{iωt}, which will be suppressed henceforth. In this case, the superposed dot denoting time derivative is replaced by i*ω* and the convolution asterisk is replaced by ordinary multiplication. Explicitly, equation (2.1) becomes
2.4
while the constitutive relations (2.2) become
2.5
If the body occupies a domain *Ω* and tractions *T*_{0}=*σ*·*n* are applied at its boundary ∂*Ω*, then the rate of working of the surface tractions and body force is
2.6
having employed the divergence theorem and the equation of motion (2.1). It follows, employing the time-reduced version of the constitutive relations (2.2), that the mean rate of working of the applied forces is
2.7
where , , etc. Also, the relevant contractions (: or ·) are left implicit. It follows, assuming *ω*>0, that positive definiteness of the quadratic forms and guarantees positive energy dissipation. This will be assumed henceforth, though relaxation of these conditions (such as ) will be discussed after the main development is complete.

Now consider the time-harmonic problem, as described above, but with more general boundary conditions. Three compatible boundary conditions are given at each point of ∂*Ω*. These could be for instance: all three components of displacement; all three components of traction; two components of traction and the orthogonal component of displacement (e.g. normal component of displacement and shear components of traction).

The problem with any such conditions can be described by the stationary principle of Willis (1981), whose time-reduced form is
2.8
the variation being taken over displacement fields *u* that satisfy whatever conditions are prescribed for displacements on ∂*Ω*. The fields *σ*_{0} and *p*_{0} are not unique but are any fields that satisfy the equation of motion (2.1), together with any conditions that are prescribed for the tractions on ∂*Ω*. To confirm this claim, note that the variation is
2.9
Requiring the volume integral to vanish for all *δu* implies the time-reduced version of the equation of motion (2.4). The surface integral picks out only those components of surface traction for which the corresponding component of displacement is not satisfied, and hence the boundary conditions on traction are enforced.

Willis (1981) also derived a dual stationary principle, whose time-reduced form is
2.10
where *u*_{0} is any displacement field that satisfies any given boundary conditions on displacement. The variation is taken over fields *σ*, *p* that satisfy the equation of motion and any given traction boundary conditions. Corresponding stationary principles for electromagnetic waves were given by Willis (1984).

## 3. Minimum variational principles for elastodynamics

First, the stationary principle (2.8) is spelled out explicitly,
3.1
It can be verified immediately that both the real and imaginary parts of equation (3.1) deliver the same Euler–Lagrange equations. Note further that equation (2.8) remains true if it is multiplied throughout by *e*^{iθ}. This is equivalent to replacing *C*, *ρ*, *σ*_{0}, *p*_{0} by *Ce*^{iθ}, *ρe*^{iθ},*σ*_{0}*e*^{iθ}, *p*_{0}*e*^{iθ} for any *θ*. Henceforth, just the imaginary part of equation (3.1) will be considered, while bearing in mind that *C*, *ρ*, *σ*_{0} and *p*_{0} can be changed as just indicated. For the validity of the reasoning to follow, *θ* should be chosen such that the imaginary parts of *Ce*^{iθ} and −*ρe*^{iθ} are positive definite. If *ρ* is real, as is typically the case, and is strictly positive definite, such a replacement allows the ensuing analysis to apply: we assume that both and are strictly positive definite, aside from in the variational principles (3.26) and (3.28).

The functional expressed in the imaginary part of equation (3.1) is saddle shaped. An alternative functional is obtained by performing a partial duality transformation, which involves and in place of and . To clarify what is going on, note first that the time-reduced constitutive relations (2.5) imply that
3.2
Now following Cherkaev & Gibiansky (1994), elementary calculation gives
3.3
where
3.4
is positive definite, as can be seen from the positivity of the associated quadratic form
3.5
The infimum is attained when satisfies (3.2)_{1}, or equivalently Re(*σ*)=Re(*Ce*). Similarly,
3.6
where
3.7
is positive definite. Attainment of this infimum is consistent with equation (3.2)_{2} or equivalently Im(*p*)=Im(i*ωρu*).

The imaginary part of the functional whose variation is indicated in equation (3.1) can now be written as
3.8
where
3.9
Furthermore,
3.10
Now,
3.11
so long as
3.12
and satisfies any given traction conditions on ∂*Ω*. Otherwise, the supremum is infinite. Therefore, the variational principle (3.9) reduces to
3.13
where is now given by equation (3.12) and the infimum is subject to any prescribed traction boundary conditions on and any prescribed displacement boundary conditions on . Since the products and do not affect the minimizer, the variational principle further simplifies to
3.14
Note that the constitutive relations (2.5) imply the new constitutive relations
3.15
The fields appearing on the left-hand side of equation (3.15) are linked by the differential constraints
3.16
while the fields on the right-hand side are linked by
3.17
Following Milton *et al.* (2009), we can rewrite equation (3.15) as a single equation , where
3.18
and *G* satisfies the differential constraints (3.16), while *F* satisfies the differential constraints (3.17) and is positive definite. If we regard equations (3.15)–(3.17) as a primal problem, then the variational principle for it is
3.19
which is equivalent to equation (3.14). To check this variational principle directly, and thus provide an alternative derivation of equation (3.14), we compute the first-order variation in the integral in equation (3.19),
3.20
This vanishes since the fields satisfy equations (3.16) and (3.17) and the stated boundary conditions. It also clearly vanishes (and produces a valid variational principle) under more general boundary conditions. For example, although difficult to realize physically, one could fix the real part of the displacement and the real part of the traction at the boundary, so that and . Then, choosing , and , the variational principle (3.14) reduces to that of Milton *et al.* (2009). Alternatively, by having no constraints on and , one obtains a variational principle in which the trial fields are not required to satisfy any boundary condition, and at the minimum, and on ∂*Ω*.

If one prefers, one can rewrite the variational principle in a form where only the surface values of and appear. Using the differential constraints on , , and and integrating by parts, equation (3.14) is equivalent to 3.21

When *f*=0, the minimum of the variational principle (3.14), or equivalently equation (3.19) is
3.22
and thus depends only on the surface displacements *u* and and surface tractions *σn* and . Thus, for any given choice of trial fields and meeting the chosen boundary conditions, which imply that
3.23
equations (3.15) and (3.22) produce the inequality
3.24
in which equation (3.23) has been used to eliminate and . Since and do not appear in this equation, we may as well choose and on ∂*Ω*, in which case we are free to choose any and . The inequality (3.24) bounds the possible surface fields *u* and *σn*, i.e. it bounds the Dirichlet to Neumann map associated with *Ω*. Alternatively if the surface fields are known from measurements, then, for each set of trial fields, the inequality provides bounds on the moduli *C* and *ρ* inside *Ω* and could potentially be used for tomography. We are assured that the inequality (3.24) will produce useful bounds. In the quasistatic limit where *ω*=0, and , the term involving 𝒫 in equation (3.24) vanishes. Then, in a periodic composite with unit cell *Ω*, equation (3.24) implies the variational inequality of Cherkaev & Gibiansky (1994) on the effective viscoelastic tensor that has provided tight bounds (Gibiansky & Lakes 1993, 1997; Gibiansky & Milton 1993; Milton & Berryman 1997; Gibiansky *et al.* 1999).

In the limit as , the last term in the variational principle (3.14) will remain finite if and only if
3.25
and in this limit, the variational principle becomes
3.26
where now *p*′′ and *u*′ are given by equations (3.12) and (3.25).

These variational principles are particularly useful if a multi-phase medium occupies *Ω* and one of the phases, occupying a subregion *Ψ* is lossless, with both *C*′′=0 and *ρ*′′=0. Then, assuming *C*′ and *ρ*′ are positive in *Ψ*, one does not have the flexibility to replace *C* and *ρ* by *Ce*^{iθ} and *ρe*^{iθ} since the imaginary parts of *Ce*^{iθ} and −*ρe*^{iθ} would lose their simultaneous positive semi-definiteness within *Ψ*. In the limit as *C*′′→0 in *Ψ*, we see from equation (3.5) that the last term in equation (3.26) will remain finite if and only if
3.27
In this limit, the variational principle becomes
3.28
in which *p*′′ and *u*′ are given by equations (3.12) and (3.25), and the infimum is over fields , such that *u*′ is an exact solution of the elastodynamic equations in *Ψ*, i.e. equations (3.17), (3.25) and (3.27).

## 4. Related elastodynamic principle of Hashin–Shtrikman type

Recall that the constitutive relations (2.5) are equivalent to those given in equation (3.15). Now introduce a ‘comparison’ medium, characterized by matrices and ,^{1} introduce ‘polarizations’ , , and , and define
4.1

These relations are consistent with the desired constitutive relations (3.15) if the polarizations are chosen so that 4.2 With this preparation, note that 4.3 Similarly, 4.4 Since equations (4.3) and (4.4) are essentially the square of the difference between the left- and right-hand sides of equation (4.2); here, ≈ means equal to a good approximation when the relations (4.2) are satisfied to a fair approximation. As discussed shortly, the ≈ sign can be replaced by an inequality when and are appropriately chosen. Equality is achieved when relations (4.2) are exactly satisfied.

The Hashin–Shtrikman variant of the minimum principle (3.13) is now 4.5 This produces an approximate solution for any given choice of the ‘polarizations’ , , and , and the exact solution if these are chosen to make the functional stationary. When and are chosen so and , then the ≈ signs in equations (4.3) and (4.4) can both be replaced by ≤ signs and one has a minimum principle that involves taking the infimum over the ‘polarizations’ , , and , as well as an infimum over the fields and . If, on the other hand, and are chosen so and , then the ≈ signs can both be replaced by ≥ signs and one obtains a principle that involves taking the supremum over the ‘polarizations’ and an infimum over the fields and . The functional in this case has only a single stationary point, a saddle, and the infimum and supremum can be taken in any order. In either case, the infimum over the fields and can be computed in terms of the relevant Green function, and one is left with a minimum principle or a maximum principle. When the body is infinite in extent, Green’s function can be explicitly computed, as shown in appendix A.

In the notation of equations (3.18) and (3.19), equation (4.5) becomes
4.6
where
4.7
and in this form, the principle is easily generalized to acoustics and electromagnetism, as will become clear in §§5 and 6. It is perhaps worth remarking that there is no need for to be restricted to have the block-diagonal structure displayed by in equation (3.18)_{1}; this, however, is not developed further at present.

Developing equation (4.6) a little further, let *F*_{0} solve the variational problem for the comparison medium, i.e. in the case that *T*=0, and write
4.8
Any admissible must satisfy homogeneous boundary conditions, and hence
4.9
It follows that equation (4.6) can be written as
4.10
The that attains the minimum must satisfy
4.11
Furthermore, is an admissible variation. Formally, writing this field as
4.12
and taking account of equation (4.11), equation (4.10) becomes
4.13
this is stationary with respect to *T* when
4.14
The operator relevant to an infinite medium is given explicitly in appendix A.

## 5. Variational principles for acoustics

For acoustics at fixed frequency in the presence of a body force *f*, the fields satisfy the differential constraints
5.1
where *P* is the pressure, *p* the momentum, *v* the velocity and *h* defined to be the divergence of the velocity. They also satisfy the constitutive relations
5.2
where *k*=1/*κ* is the compressibility, the inverse of the bulk modulus *κ*, which is complex with *k*′′<0, and *r*=*ρ*^{−1} is the inverse density (which we allow to be matrix valued and complex with *r*′′>0). Introducing the positive definite matrices
5.3
the constitutive relations take the new form
5.4
which can be expressed as the single equation , with
5.5
where the fields appearing in *G* are linked by the differential constraints
5.6
while the fields appearing in *F* are linked by
5.7
(The second part in equation (5.4) was written with the real components of the fields on the left, rather than on the right, to ensure that the differential constraints couple fields on the same sides of equation (5.4).) Based on equation (3.19), we assert that the form of the variational principle for acoustics is
5.8
in which *h*′′ and *p*′′ are given by equation (5.7), and *h*′_{0}, *P*′′_{0}, *v*′_{0} and *p*′_{0} satisfy the differential constraints (5.6). To check this we see that at the solution of (5.4)–(5.7), the first-order variation
5.9
only depends on boundary terms. (Note the need for the factor *ω*^{2} in equation (5.9), which explains the need for the additional factors of *ω* in the last two components of *G* and *F* in equation (5.5).)

The condition that the boundary integral vanishes for all *δP*′ and *δv*′′ that are permitted according to which components of *P*′ and *v*′′*n* are prescribed on ∂*Ω*, forces the complementary components of (*v*′_{0}−*v*′)*n* and *P*′′_{0}−*P*′′ to be zero, thus fixing these components of *v*′*n* and *P*′′. In the particular case when *f*=0 and all components of *P*′ and *v*′′*n* are prescribed at the boundary, then by choosing *P*′′_{0}=0, *p*′_{0}=0 and *v*′_{0}=0, these variational principles reduce to those of Milton *et al.* (2009).

If one prefers, one can rewrite equation (5.8) in a form where only the surface values of *P*′′_{0} and *v*′_{0}*n* appear. Using integration by parts, the variational principle is equivalent to
5.10

When *f*=0, the minimum value of equation (5.8) or equation (5.10) only depends on the boundary fields, and is
5.11
Thus, knowledge of this boundary data will, for any given choice of trial fields and , give constraints
5.12
on the possible values of the moduli *κ* and *ρ* inside *Ω*, and thus may be useful for tomography.

In the limit as *r*′′→0, i.e. as *ρ*′′→0, the last term in equation (5.8) will be minimized when
5.13
and is infinite otherwise. In this limit, the variational principle reduces to
5.14
where *h*′′ and *p*′′ are given by equation (5.7), while *v*′′ is given by equation (5.13).

## 6. Variational principles for electrodynamics

For electromagnetism, it is more conventional to have fields that are multiplied by *e*^{−iωt}, rather than *e*^{iωt}, so that the permittivity tensor *ε* and permeability tensor *μ* have positive definite imaginary parts. We now adopt this convention. Then, the fields in Maxwell’s equations at fixed frequency *ω* satisfy the differential constraints
6.1
and the constitutive equations
6.2
where *m*=*μ*^{−1}. Introducing the positive definite matrices
6.3
the constitutive relations take the new form
6.4
which can be expressed as the single equation , with
6.5
where the fields appearing in *G* are linked by the differential constraints
6.6
while the fields appearing in *F* are linked by
6.7
Based on equation (3.19), we state that the form of the variational principle for electromagnetism is
6.8
where *D*′ and *B*′′ are given by equation (6.7), while *D*′′_{0}, *E*′′_{0}, *H*′_{0} and *B*′_{0} satisfy the differential constraints (6.6). To check this, we use the identities
6.9
which hold for any fields satisfying the differential constraints (6.6) and (6.7). Then, we see that at the solution of (6.4)–(6.7), the first-order variation can be expressed in terms of boundary fields
6.10
Note that this surface integral only depends on the tangential components of the fields *H*′−*H*′_{0}, *δE*′, *E*′′_{0}−*E*′′ and *δH*′′. The condition that the boundary integral vanishes for all *δE*′ and *δH*′′ that are permitted according to which tangential components of *E*′ and *H*′′ are prescribed, forces the complementary tangential components of *H*′−*H*′_{0} and *E*′′_{0}−*E*′′ to be zero, thus fixing these tangential components of *H*′ and *E*′′. Again, these variational principles reduce to those of Milton *et al.* (2009) when the tangential components of *E*′ and *H*′′ are fully prescribed.

If one prefers, one can rewrite equation (6.8) in a form where only the surface values of *H*′_{0} and *E*′′_{0} appear. Using integration by parts, the variational principle is equivalent to
6.11

When *j*=0, the minimum value of equation (6.8) or equation (6.11) is
6.12
Consequently, for any given choice of trial fields and , we have the inequality
6.13
which could be useful for tomography if we have measurements of the boundary values of the fields.

In the limit as *μ*′′→0, i.e. as *m*′′→0, the last term in equation (6.8) will be minimized when
6.14
and is infinite otherwise. In this limit, the variational principle reduces to
6.15
where *D*′ and *B*′′ are given by equation (6.7), while *H*′′ is given by equation (6.14). When, in addition, *ε*′′→0 in some subregion *Ψ* of *Ω*, there is a variational principle analogous to equation (3.28),
6.16
where *D*′, *B*′′ and *H*′′ are given by equations (6.7) and (6.14), while the infimum is over fields *E*′, which are an exact solution of Maxwell’s equations inside *Ψ*, i.e. satisfying *D*′=*ε*′*E*′ within *Ψ* as well as equations (6.7) and (6.14).

## Acknowledgements

The authors are grateful to the referees for their many helpful comments and are thankful for support from the National Science Foundation through grant DMS-0707978.

## Appendix A. Green’s function for infinite comparison medium

The Hashin–Shtrikman variational functional is rendered stationary (with respect to the fields , , , ) when is the strain field associated with the displacement and the equation of motion A1 is satisfied, together with the relations (4.1).

Motivated by the form of relations (3.2), define
A2
in which *D*_{2}, *D*_{3}, *Q*_{2} and *Q*_{3} are symmetric real matrices and, in terms of these and the real matrices *D*_{1} and *Q*_{1},
A3
The relations (4.1) may then be given in the equivalent form
A4

The quadratic forms associated with and are positive definite if *D*_{2},*D*_{3},−*Q*_{2} and −*Q*_{3} are positive definite, as will henceforth be assumed.

The equation of motion (A1) now gives A5 where A6

The required Green function is that associated with equation (A5). It can be constructed by a method introduced for the ‘ordinary’ elastodynamic Green function by Willis (1980). First, select , with constant, and note the plane-wave decomposition of the three-dimensional delta function
A7
This motivates first solving equation (A6), with replaced by , where |*ξ*|=1. The solution, *U*^{ξ} say, depends on *x* only in the combination *s*=*ξ* · *x*; thus, it satisfies the system of ODEs
A8
where
A9
The system can be diagonalized by reference to the eigenvalue problem
A10
Although the matrices are real and symmetric, they are not positive definite, and the eigenvalues may be real or they may occur in complex conjugate pairs. Suppose first that equation (A10) is satisfied for some real and (hence real) *U*_{N}, having the form shown in equation (A6)_{1}: say , where and are 3-vectors. The 3-vector can be eliminated from equation (A10) to give
A11
and hence
A12
Recalling the definiteness of *D*_{2},*D*_{3},−*Q*_{2} and −*Q*_{3}, it follows that the expression on the left-hand side of equation (A12) cannot be zero if is positive. Hence, any real eigenvalue must be negative. Now write
A13
Since cannot be zero, it is possible to ensure that for all *N*. If, furthermore, is complex, then and for some *M*, and correspondingly . The eigenvectors are normalized so that
A14
Setting
A15
it follows that
A16
The solution that decays to zero as is
A17
Having now determined *U*^{ξ}(*s*), the solution *U*(*x*) of equation (A5), with , follows by making the superposition implied by equation (A7),
A18
where
A19
If, for any *N*, the eigenvalue and eigenvector are real, then and the contribution to the sum is real. If the eigenvalue is complex, it follows from the properties established above that some other *M* contributes the complex conjugate expression. Thus, is real.

Now for some general , and in particular the given by equation (A6), the solution is A20 which gives us and and hence and , and from equation (A4), we get A21

The functions , , and thus obtained in terms of the polarization fields are then substituted back in equation (4.5). The resulting functional is a minimization (or maximization) principle for the polarization fields when and are chosen so that both and are positive definite (respectively, negative definite).

This appendix is concluded by recording the explicit form of the operator introduced in equation (4.12). For this purpose, it is helpful to define
A22
so that the strain *e* associated with displacement *u* is , and also to define
A23
but with the partial derivatives acting backwards with respect to the immediately preceding variable, so that, for example,
A24
We also write
A25
In terms of the representation (4.12) employing the polarization *T* as defined in equation (4.7), it follows from equations (A5), (A20) and (A21) that
A26
The symmetry of is consistent with its association with the variational formulation for the comparison medium.

## Footnotes

- © 2010 The Royal Society