Large-Scale Structures in Acoustics and Electromagnetics: Proceedings of a Symposium (1996)

Adrianus T. de Hoop

Delft University of Technology Delft, The Netherlands

Maarten V. de Hoop

Schlumberger Cambridge Research Cambridge, England

Acoustic, elastodynamic, and electromagnetic wavefields share a number of common features in their mathematical description. Their local, pointwise behavior in space-time is governed by a hyperbolic system of first-order partial differential equations that are representative for the physical phenomena involved on a local scale. When supplemented with initial conditions that relate the wave solution to its excitation mechanism and boundary conditions across interfaces where the coefficients in the system jump by finite amounts, the problem has a unique solution. The computational handling of the problem, however, often starts from a "weak" formulation, where the pointwise, or "strong," satisfaction of the equality signs in the equations is replaced by requirements on the equality of integrated or weighted

versions of the differential equations. The resulting expressions have a counterpart in physics that is found in the pertaining reciprocity theorems of the Rayleigh (acoustic waves in fluids), Betti-Rayleigh (elastic waves in solids), or H.A. Lorentz (electromagnetic waves) types. This observation has led to the approach presented in the present paper, where the relevant reciprocity theorems are taken as points of departure. Through them, a computational scheme is conceptually taken to describe the interaction between the actual wavefield state to be computed and a suitably chosen "computational state" that is representative for the method at hand, just as in physics the reciprocity theorems describe the interaction between observing state and observed state, or quantify the reciprocity between transmitting and receiving properties of any device or system (transducer in acoustics and elastodynamics, antenna in electromagnetics, electromagnetic compatibility of interfering electromagnetic systems or devices). It is also believed that through this point of view one is guided to developing computational algorithms for each of the three types of wave fields in a manner that expresses the structures common to all of them. Background literature on reciprocity can be found in some papers by A.T. de Hoop (1987, 1988, 1989, 1990, 1991, 1992) and in a forthcoming book (de Hoop, 1995).

We consider linear acoustic, elastic, or electromagnetic wave motion in some subdomain of three-dimensional Euclidean space . The configuration in which the wave motion is considered to be present is assumed to be time invariant and linear in its physical behavior. The wave quantities involved are found to satisfy certain reciprocity properties which will be taken as the point of departure for our further considerations. Now, for the indicated type of configuration, there prove to be two kinds of reciprocity theorem: one of the time-convolution type, the other of the time-correlation type. Several operations on the wave quantities will occur throughout the paper. First, we shall introduce their notation.

Cartesian coordinates x={x₁,x₂,x₃} are used to specify position; t is the time coordinate. Differentiation with respect to x_p is denoted by is a reserved symbol for differentiation with respect to t. The subscript notation for the vectorial and tensorial quantities occurring in the wave motion will be used whenever appropriate; the subscripts are to be assigned the values 1, 2, and 3.

The characteristic function of the domain is denoted by χ_D and is given by

where is the boundary of , and is the complement of in

Let F = F(x,t) denote any space-time function. Then, the time reversal operator T is defined by

It has the property

Let Q(x,t) denote another space-time function, then the time convolution C_t(F,Q) of F and Q is defined as

It has the properties

The time correlation R_t (F, Q) of F and Q is defined as

It has the properties

Let the one-dimensional field matrix F_p = F_p(x,t) of the wave motion be composed of the components of the two wavefield quantities whose inner product represents the area density of power flow (Poynting vector). Then, F_p satisfies a system of linear, first-order, partial differential equations of the general form

where uppercase Latin subscripts are used to denote the pertaining matrix elements and the summation convention for repeated subscripts applies. In (3.12), D_I,P is a symmetrical, block off-diagonal spatialdifferentiation operator matrix that contains the operator in a homogeneous linear fashion that is specific for each type of wave motion under consideration, M_I,P = M_I,P (x) is the medium matrix that is representative for the physical properties of the (arbitrarily inhomogeneous, anisotropic) medium in which the waves propagate, and Q_I = Q_I (x,t) is the volume source density matrix that is representative for the action of the volume sources that generate the wavefield.

The medium parameters are assumed to be piecewise continuous functions of position. Across a surface of discontinuity in medium properties, the parameters may jump by finite amounts. On the assumption that the interface is passive (i.e., free from surface sources) and that the wavefield quantities

must remain bounded on either side of the interface, the wavefield must satisfy the boundary condition of the continuity type

where N_I,P is the unit normal operator at the interface that arises from replacing in D_I,P by n_p where n_p is the unit vector along the normal to the interface.

For acoustic waves in fluids,

where p = acoustic pressure and v_r = particle velocity, and

where q = volume source density of injection rate and f_k = volume source density of force. For elastic waves in solids,

where v_r = particle velocity and τ_p,q = dynamic stress, and

where f_k = volume source density of force and h_i,j = volume source density of deformation rate. For electromagnetic waves,

where E_r = electric field strength and H_p = magnetic field strength, and

where J_k = volume source density of electric current and K_j = volume source density of magnetic current. The structures of D_I,P and M_I,P for the three types of wave fields are given in Appendix 3A.

In the reciprocity theorems to be discussed below, two diagonal matrices and occur that concatenate, out of the wavefields pertaining to two admissible states, their interaction. For acoustic waves in fluids the diagonal matrix is given by

for elastic waves in solids by

and for electromagnetic waves by

The diagonal matrix is just the unit matrix:

For the reciprocity theorem of the time-convolution type to hold, a necessary and sufficient condition proves to be

This condition requires that the block-diagonal part of D_I,P be anti-symmetric and that its block off-diagonal part be symmetric. For the reciprocity theorem of the time-correlation type to hold, a necessary and sufficient condition proves to be

This condition requires that D_I,P be symmetric. The two conditions are independent, but if they are satisfied simultaneously, D_I,P is a symmetric, block off-diagonal matrix operator. For the three types of wave motion considered in this paper, this is indeed the case. It is therefore conjectured that the indicated structure of the spatial differential matrix operator could prove to be fundamental in order that a system of first-order partial differential equations be representative for a physical wave motion.

It is noted that the medium matrix M_I,P is not subjected to any restriction of this kind.

In view of the linearity of the wave motion, the principle of superposition ensures that the wave field F_p that is generated by the volume source distribution Q_I can be written as the superposition of point-source contributions through the use of a Green's tensor. The latter is a solution of the system of differential equations

where is the unit matrix and δ(x-x',t-t') is four-dimensional Dirac delta distribution operative at {x,t}= {x',t'}. In view of the time invariance of the medium, the Green's tensor depends on t and t' only through the difference t - t', i.e., G_P,I' = G_P,I'(x,x',t,t') = G_P,I'(x,x',t - t'). The Green's tensor plays an important role in the embedding formulations of the wave field problem.

In the wavefield reciprocity theorems certain interaction quantities are considered that are representative for the interaction between two admissible states of the pertaining wavefield in a given (proper or improper) subdomain D of . Each of the two states applies to its own medium and has its own volume source distribution. Let the superscripts A and Z indicate the two states, then the wavefields in the two states are related to their respective sources via

Further, for each of the two states the boundary condition of the continuity type

holds.

The local interaction quantity to be considered in the reciprocity theorem of the time-convolution type is , where the property of (3.14) has been used. With the aid of (3.17) and (3.18) this expression is rewritten as

Equation (3.21) is the local form of the reciprocity theorem of the time-convolution type. The global form, for the domain , of this theorem follows upon integrating (3.21) over the domain and applying Gauss' integral theorem to the first term on the left-handed side over each subdomain of where the field quantities are continuously differentiable. Adding the contributions from these subdomains, the contributions from sourcefree interfaces of discontinuity in medium properties in the interior of cancel in view of the boundary conditions given in (3.19) and (3.20) and only a surface integral over the boundary of remains. The result is

Equation (3.22) is the global form, for the domain , of the reciprocity theorem of the time-convolution type.

The terms in (3.21) and (3.22) containing the medium matrices define the contrast-in-medium contributions to the time-convolution interaction of the two states. They vanish at those positions where . If this condition holds, the media in the two states are denoted as each other's adjoints. If the condition holds for one and the same medium, such a medium is denoted as self-adjoint. An isotropic medium is always self-adjoint. The terms containing the volume source densities yield the contribution from the volume sources to the interaction of the two states. They vanish at sourcefree positions.

In a number of applications (3.22) will be applied to the entire . Then, outside some sphere S(O, Δ₀) with radius Δ₀ and center at the origin O of the chosen reference frame, the media in the two states will be assumed to be the same and homogenous as well as isotropic. For such a medium, the tensor Green's function is known analytically and in particular the causal and anti-causal source-type integral representations are known analytically. For the application of (3.22) to the entire , the theorem will be first applied to a sphere S(O, Δ) of radius Δ and center at the origin O of the chosen reference frame and the limit will be taken. If, now, in both states the wavefields are causally related to the action of their volume source distributions (assumed to have bounded supports), the integral over S(O, Δ) vanishes as . However, if one of the two states is causally related to the action of its volume sources and the other anti-causally, the integral over S(O,Δ) does not vanish as , but has a constant value for sufficiently large values of Δ.

The local interaction quantity to be considered in the reciprocity theorem of the time-correlation type is , where the property of (3.15) has been used. With the aid of (3.11), (3.17), and (3.18) this expression is rewritten as

Equation (3.23) is the local form of the reciprocity theorem of the time-correlation type. The global form, for the domain , of this theorem follows upon integrating (3.23) over the domain and applying Gauss' integral theorem to the first term on the left-hand side over each subdomain of where the field quantities are continuously differentiable. Adding the contributions from these subdomains, the contributions from interfaces of discontinuity in medium properties in the interior of cancel in view of the boundary conditions given in (3.19) and (3.20) and only a surface integral over the boundary of remains. The result is

Equation (3.24) is the global form, for the domain , of the reciprocity theorem of the time-correlation type.

The terms in (3.23) and (3.24) containing the medium matrices define the contrast-in-medium contributions to the time-correlation interaction of the two states. They vanish at those positions where . If this condition holds, the media in the two states are denoted as each other's time reverse adjoints. (The ''time-reverse'' is reminiscent of the fact that "adjoint" applies to the reciprocity theorem of the time-convolution type and that correlation can be considered as a compound operation consisting of convolution and time reversal.) If the condition holds for one and the same medium, such a medium is denoted as time-reverse self-adjoint. For an isotropic medium, the medium matrix is diagonal; an instantaneously reacting isotropic medium is therefore always time-reverse self-adjoint. The terms containing the volume source densities yield the contribution from the volume sources to the interaction of the two states. They vanish at source-free positions.

In a number of applications (3.24) will be applied to the entire . Then, outside some sphere S(O,Δ₀) with radius Δ₀ and center at the origin O of the chosen reference frame, the media in the two states will be assumed to be the same and homogenous as well as isotropic. For such a medium, the tensor Green's function is known analytically and in particular the causal and anti-causal source-type integral representations are known analytically. For the application of (3.24) to the entire , the theorem will be first applied to a sphere S(O, Δ) of radius Δ and center at the origin O of the chosen

reference frame and the limit will be taken. If, now, in State A the wave field is causally related to the action of its volume source distributions, and in State Z the wavefield is anti-causally related to the action of its volume source distributions (the volume source distributions being assumed to have bounded supports), the integral over S(O,Δ) vanishes as . However, if both states are causally related to the action of their volume sources, the integral over S(O,Δ) does not vanish as , but has a constant value for sufficiently large values of Δ.

For the choice State A = State Z and zero correlation time shift (i.e., t = 0), (3.23) reduces to the local energy balance for the wavefield and (3.24) to the global energy balance for the domain , provided that M_Q.P = M_P,Q. This implies that for the energy considerations pertaining to a physical wavefield to hold, the medium matrix must be symmetric. In that case, also the quantity (1/2)M_P.QF_PF_Q (whose time derivative occurs in equations (3.23) and (3.24)) should represent the volume density of stored energy. For the latter, the symmetric medium matrix should, in addition, on physical grounds be positive definite.

Equation (3.22) leads to a reciprocity property of the Green's tensor. Let be the causal wavefield in Medium A generated by the point source operative at {x,t} = {x',0}. Then (cf. (3.16)) . Let, similarly, be the causal wavefield in Medium Z generated by the point source operative at {x,t} = {x",0}. Then . Take the media in the two states as each other's adjoints, i.e., , and apply (3.22) to the entire . In this application, the contrast-in-media term and the contribution from the "sphere at infinity" vanish. The result is

Since (3.25) has to hold for arbitrary values of and we end up with

Equation (3.26) is the reciprocity relation for the causal Green's tensor.

On many occasions the wavefield computation in an entire configuration is beyond the capabilities because of the storage capacity required and the computation time involved. In that case, it is standard practice to select a target region of bounded support in which a detailed computation is to be

carried out, while the medium in the remaining part of the configuration (the embedding) is chosen to be so simple that the wave motion in it can be determined with the aid of analytical methods. In particular, this applies to scattering problems and to geophysical modeling, where the support of the model configuration is often taken to be the entire . Examples of such simple embeddings are a homogeneous, isotropic medium (chosen for most scattering configurations) and a medium consisting of parallel layers of homogeneous, isotropic material (chosen in most geophysical applications). In these cases, time Laplace and spatial Fourier transform techniques provide the analytical tools to determine the wave motion or, in fact, the relevant Green's tensor. Once the embedding has been chosen, the problem of computing the wavefield in can be formulated as a contrast problem . For this, we proceed as follows.

The State A is introduced consisting of the Actual wavefield , the actual sources that excite it, and the actual medium in which the propagation takes place. Next, we introduce a State B consisting of the wavefield that the actual sources would generate in the medium of the emBedding. Denoting the Green's tensor of the embedding by the latter wavefield is expressible as

where is the support of the sources that generate the wavefield in the actual configuration. From the corresponding wavefield equations it then follows that

This equation can be rewritten in two ways as a wavefield equation for the Contrast state, to be denoted by the superscript C, in which the contrast wave field is

In one of them, the medium properties in the wave operator on the left-hand side are taken to be the ones of the actual medium; this is typically done in the combination of the embedding technique with finite-element or finite-difference modeling. In the other, the medium properties in the wave operator on the left-hand side are taken to be the ones of the embedding; this is typically done in the integral-equation or method-of-moments modeling. Both ways lead to a contrast-source formulation. The expressions for the two cases are given below.

For the contrast formulation for finite-element or finite-difference modeling, (3.28) is, in combination with (3.29), rewritten as

with

and

Note that in this contrast formulation, the contrast source density is known.

For the contrast formulation for integral-equation or methods-of-moments modeling, (3.28) is, in combination with (3.29), rewritten as

with

and

Note that in this contrast formulation, the contrast source density is unknown, since is unknown.

In the two sections following, it will be indicated how these different states are used in the reciprocity theorems of the previous section to lead to computational schemes for the evaluation of the wavefields.

In finite-element/finite-difference modeling over an entire configuration occupying the bounded domain , the wavefield to be computed is the total wavefield. The latter is approximated by an expansion of the type

where is an appropriate sequence of known, linearly independent expansion functions with as their supports, and is the sequence of expansion coefficients to be computed. In typical finite-element/finite-difference modeling the support of each expansion function is an elementary subdomain of the (discretized) version of (usually a simplex or a complex). Further, boundary conditions as needed for the uniqueness of the solution in are prescribed on . Next, a sequence of "computational" states, denoted by the superscript Z, is selected, for which

where the right-hand side is a sequence of known, linearly independent weighting functions with as their supports. Finally, we take

and hence

Application of the earlier reciprocity theorems to the State A and the sequence of States Z leads to a system of linear algebraic equations in the expansion coefficients.

In case an embedding procedure is applied, the wavefield to be computed is the contrast wavefield. The latter is approximated by an expansion of the type

where is an appropriate sequence of known, linearly independent expansion functions with as their supports and is now the sequence of expansion coefficients to be computed. In typical finite-element/finite-difference modeling the support of each expansion function is an elementary subdomain of the (discretized) version of (usually a simplex or a complex). Further, "absorbing boundary conditions" as needed for the uniqueness of the solution in are prescribed on . These should model the radiation of the contrast wavefield into the passive embedding. Next, a sequence of "computational" states, denoted by the superscript Z, is selected, for which

where the right-hand side is a sequence of known, linearly independent weighting functions with as their supports. Finally, we take

and hence

Application of the reciprocity theorems to the State A and the sequence of States Z again leads to a system of linear algebraic equations in the expansion coefficients.

In finite-element/finite-difference modeling, the expansion and weighting functions are standardly taken to be polynomials in the time variable and the spatial variables. Their vector and tensor components in space can be organized such that the continuity conditions across an interface between two different media are taken into account automatically, while leaving those components that are not necessarily continuous free to jump by finite amounts. Such a procedure can be carried out consistently if the simplex is taken as the elementary subdomain of the discretized configuration and a consistent linear approximation within each simplex is used. Thus, the notions of "face element" and "edge element" for arbitrary vectors and tensors have been introduced. For literature on the subject, see Mur and de Hoop (1985) and Mur (1990, 1991, 1993) for the application to electromagnetic fields and Stam and de Hoop (1988, 1989, 1990) for the application to elastodynamic wavefields.

The integral-equation or method-of-moments modeling is invariably based on an embedding procedure. As a consequence of this, the wavefield to be computed is the contrast wavefield. The latter is approximated by

where is an appropriate sequence of known, linearly independent expansion functions with as their supports, and is the sequence of expansion coefficients to be computed. The contrast source density is written as

in which the first term on the right-hand side is known and the second term on the fight-hand side is unknown. Next, a "computational" state, denoted by the superscript Z, is selected, for which

and

Substitution in the earlier reciprocity theorems then leads to a system of linear, algebraic equations in the expansions coefficients.

Although the real, physical wave phenomena take place in space-time, it can under certain circumstances be advantageous to parametrize the problem in the coordinates in which shift invariance in the configuration occurs. Since we have assumed that our configurations are, apart from linear, time-invariant in their physical behavior, such a procedure certainly applies to the time coordinate. Moreover, in this coordinate the principle of causality applies. In view of these two aspects, the time Laplace transformation performs the appropriate parametrization in the time coordinate. For any causal, bounded function Q_I = Q_I(x,t) with temporal support this transformation is

Here, s is the time Laplace transform parameter or complex frequency. The time Laplace transformation has the following properties:

In view of Lerch's theorem (Widder, 1946), the correspondence between and Q_I(x,t) for t>t₀ is unique. Using these properties, the space-time wave motion can be recovered after having solved a sequence of space problems with appropriate values of the time Laplace transform parameter. For recent results in this direction, see Lee et al. (1994).

In the preceding two sections it has been indicated how finite-difference/finite-element methods and integral-equation/method-of-moments methods for the computation of wavefields can be envisaged to arise from the time-convolution and time-correlation type reciprocity theorems pertaining to these wavefields. This does not mean that all possibilities in this respect have found application as yet. Apart from the different choices that can still be made in the selection of the sequences of expansion and weighting functions, it also happens that, for example, the application of the reciprocity theorem of the time-correlation type to the integral-equation modeling of forward wave scatttering problems has, as far as the present authors are aware, not been pursued yet. This is the more remarkable since this theorem finds prime application in the modeling of inverse scattering problems and has been extensively used in this realm. Whether or not the missing applications in a total matrix of possibilities might lead to better algorithms remains to be investigated.

The first author (A.T. de Hoop) acknowledges with gratitude the financial support from the Stichting Fund for Science, Technology, and Research (a companion organization to the Schlumberger Foundation in the United States) for carrying out the research presented in this paper.

de Hoop, A.T., 1987, "Time-domain reciprocity theorems for electromagnetic fields in dispersive media,"Radio Science22(7), 1171-1178.

de Hoop, A.T., 1988, "Time-domain reciprocity theorems for acoustic wave fields in fluids with relaxation,"The Journal of the AcousticalSociety of America84(5), 1877-1882.

de Hoop, A.T., 1989, "Reciprocity theorems for acoustic wave fields in fluid/solid configurations,"The Journal of the Acoustical Societyof America87(5), 1932-1937.

de Hoop, A.T., 1990, "Reciprocity, discretization, and the numerical solution of elastodynamic propagation and scattering problems." In: Elastic Waves and Ultrasonic Nondestructive Evaluation, S.K. Datta, J.D. Achenbach, and Y.S. Rajapakse (eds.), Amsterdam: Elsevier Science Publishers, 87-92.

de Hoop, A.T, 1991, "Reciprocity, discretization, and the numerical solution of direct and inverse electromagnetic radiation and scattering problems,"Proceedings of the IEEE79(10), 1421-1430.

de Hoop, A.T., 1992, "Reciprocity, causality, and Huygens' principle in electromagnetic wave theory." In: Huygens' Principle 1690-1990:Theory and Applications, H. Blok, H.A. Ferwerda, and H.K. Kulken (eds.), Amsterdam: Elsevier Science Publishers, 171-192.

de Hoop, A.T., 1995, Handbook of Radiation and Scattering of Waves , London: Academic Press.

Lee, K.H., G.Q. Xie, T.M. Habashy, and C. Torres-Verdin, 1994, "Wavefield transform of electromagnetic fields,"Society of Exploration GeophysicistsInternational Exposition and 64th Annual Meeting, Los Angeles, October 23-28, 1994, Expanded Abstracts, pp. 633-635.

Mur, G., 1990, "A mixed finite-element method for computing three-dimensional time-domain electromagnetic fields in strongly inhomogeneous media,"IEEE Transactions on MagneticsMAG-26(2), 674-677.

Mur, G., 1991, "Finite-element modeling of three-dimensional electromagnetic fields in inhomogeneous media,"Radio Science26(1), 275-280.

Mur, G., 1993, "The finite-element modeling of three-dimensional electromagnetic fields using edge and nodal elements,"IEEE Transactionson Antennas and PropagationAP-41(7), 948-953.

Mur, G., and A.T. de Hoop, 1985, "A finite-element method for computing three-dimensional electromagnetic fields in inhomogeneous media,"IEEE Transactions on MagneticsMAG-21(6), 2188-2191.

Stam, H.J., and A.T. de Hoop, 1988, "Time-domain reciprocity theorems for elastodynamic wave fields in solids with relaxation and their application to inverse problems,"Wave Motion10, 479-489.

Stam, H.J., and A.T. de Hoop, 1989, "A space-time finite-element method for the computation of three-dimensional elastodynamic wave fields (theory)." In: Elastic Wave Propagation, M.F. McCarthy and M.A. Hayes (eds.), Amsterdam: Elsevier Science Publishers, 483-488.

Stam, H.J., and A.T. de Hoop, 1990, "Theoretical considerations on a finite-element method for the computation of three-dimensional space-time elastodynamic wave fields,"Wave Motion12, 67-80.

Widder, D.V., 1946, The Laplace Transform, Princeton, N.J.: Princeton University Press, 63-65.

In this appendix the structures of the spatial differential operator and the medium matrix in the system of (3.12) for acoustic waves in fluids, elastic waves in solids, and electromagnetic waves are given.

For acoustic waves in fluids, the spatial differential operator in the system of (3.12) has the following form:

The medium matrix is given by

where k is the compressibility and ρ_k,r is the volume density of (inertial) mass.

For elastic waves in solids, the spatial differential operator in the system of (3.12) has the following form:

in which

and

The medium matrix is given by

where ρ_k,r is the volume density of (inertial) mass and S_t,j,p,q is the compliance.

For electromagnetic waves the spatial differential operator in the system of (3.12) is given by

The medium matrix is found to be

where ε_k,r is the permittivity and μ_J,P is the permeability.