Unified matrix-vector wave equation, reciprocity and representations
Kees Wapenaar
physics.app-ph




c
C
u
U





∇
_1
_2
k

m
p

q
d

q_1
q_2



x
x_H

S
T
V


N
K
G
J
O
I

  






   
   

ŭ
k̨^̨-̨1̨
ç
ϑ
X̋
Y
c'
a
ḇ
j
h^b
h^b
q^f
χ






ıi
i
j




𝔸
𝔻
∂
_0∪_1
1/2
1/4



E_33”
E_33'

_3”
_3'

_1”


A JKR solution for a ball-in-socket contact geometry as a bi-stable adhesive system
    M. Ciavarella
    Received ***; accepted ***
===================================================================================




The matrix-vector wave equation is a compact first-order differential equation. It was originally used for the analysis of elastodynamic plane 
waves in laterally invariant media. It has been extended by various authors for laterally varying media.
Other authors derived a similar formalism for other wave phenomena. 
This paper starts with a unified formulation of the matrix-vector wave equation for 3D inhomogeneous, dissipative media. 
The wave vector, source vector and operator matrix are specified in the appendices for 
acoustic, quantum mechanical, electromagnetic, elastodynamic, poroelastodynamic, piezoelectric and seismoelectric waves. 
It is shown that the operator matrix obeys unified symmetry relations for all these wave phenomena. 
Next, unified matrix-vector reciprocity theorems of the convolution and correlation type are derived, utilising the symmetry properties of the operator matrix. 
These theorems formulate mathematical relations between two wave states in the same spatial domain.
A unified wave field representation is obtained by replacing one of the states in the convolution-type reciprocity theorem by a Green's state.  
By replacing both states in the correlation-type reciprocity theorem by Green's states, a unified representation of the homogeneous Green's matrix is obtained.
Applications of the unified reciprocity theorems and representations for forward and inverse wave problems are briefly indicated.





Electromagnetic theory, Theoretical seismology,  Wave propagation





§ INTRODUCTION

The basic equations for wave propagation in an inhomogeneous medium can be organised in a compact matrix-vector wave equation.
This equation expresses the vertical derivative of a wave vector in terms of an operator matrix acting on this wave vector.
This specific form of the wave equation is useful, for example, to evaluate wave problems in media of which the medium parameters vary more 
rapidly in the vertical direction than in the lateral directions. 
It is also particularly useful for situations in which the vertical direction is the preferred direction of wave propagation.
However, the theoretical treatment of the matrix-vector wave equation in this paper is not limited to these special situations.







The  matrix-vector wave equation finds its roots in early work on the analysis of plane waves in laterally invariant media.
<cit.> introduced a matrix formalism for the analysis of elastodynamic plane waves propagating through a stratified solid medium.
<cit.> used the same formalism to analyse the dispersion of surface waves in layered media.
<cit.> used similar concepts to derive long-wave effective anisotropic parameters for stratified media. This approach has become known as Backus averaging <cit.>.
<cit.> used the matrix-vector wave equation to derive so-called propagator matrices for elastodynamic wave problems in stratified media.
<cit.> extended the formalism for arbitrary anisotropic inhomogeneous media and used it for the study of surface waves in  laterally varying layered media.
<cit.>, <cit.>, <cit.> and <cit.> used the matrix-vector wave equation to derive symmetry properties 
of reflection and transmission responses of  laterally invariant media.
<cit.>, <cit.>, <cit.> and <cit.> exploited the symmetry properties of the matrix-vector wave equation to 
derive so-called propagation invariants for laterally varying layered media and used this for modelling of reflection and transmission responses of such media.
Using the same symmetry properties, <cit.> and <cit.> derived  reciprocity theorems and representations for the acoustic wave vector. 

The matrix-vector wave equation has been used by many authors as the starting point for decomposition into 
coupled wave equations for downgoing and upgoing waves, 
for example for modelling in horizontally layered media  <cit.>, 
wide-angle propagation in laterally variant media <cit.>, 

reciprocity theorems for coupled down- and upward propagating waves <cit.>, 
generalised Bremmer series representations for reflection data <cit.>
and seismic interferometry <cit.>.

This paper discusses the matrix-vector wave equation and its symmetry properties for a range of wave phenomena in a unified way (section <ref>
and Appendices <ref> to <ref>).
The treatment builds on earlier systematic treatments of different wave phenomena by <cit.>, <cit.>, <cit.>, <cit.>,
<cit.>, <cit.>, <cit.>, <cit.> and <cit.>.
The matrix-vector wave equation forms the basis for the derivation of unified
matrix-vector reciprocity theorems (section <ref>)  and representations (section <ref>), 
analogous to those for the acoustic wave vector <cit.>.
 




§ THE UNIFIED MATRIX-VECTOR WAVE EQUATION AND ITS SYMMETRY PROPERTIES




 §.§ The matrix-vector wave equation





The unified matrix-vector wave equation has the form


    ∂_3 =  +.


Here  is the wave field vector,   the source vector and   the operator matrix. 
All quantities are defined in the space-frequency domain, hence =(,ω), etc., where  denotes the
Cartesian coordinate vector (x_1,x_2,x_3) and ω the angular frequency. The positive x_3-axis is pointing downward.
Operator ∂_3 stands for the spatial differential operator ∂/∂ x_3.
The vectors and matrix in equation (<ref>) are partitioned as follows


    =[ ;  ],  =[ _1; _2 ],  = [ _11 _12; _21 _22 ],


hence,


    ∂_3   =   _11+_12+_1,
    ∂_3   =   _21+_22+_2.


The vectors  and  are specified in rows <ref> to <ref>
of Table <ref> for the different wave phenomena considered in this paper. 
The wavefield quantities making up these vectors are defined in Appendices <ref> to <ref>. 
For acoustic and quantum mechanical waves  (rows <ref>  and <ref>),   and  are scalars. 
For electromagnetic, elastodynamic and poroelastodynamic waves (rows <ref>  to <ref>) they are 
2× 1, 3× 1 and 4× 1 vectors, respectively
(superscripts b, f and s in row <ref> stand for bulk, fluid and solid, respectively).
Rows <ref>  and <ref> represent coupled electromagnetic and (poro)elastodynamic waves. 

For piezoelectric waves (row <ref>), constitutive equations (<ref>) and (<ref>) account for the coupling.
For this situation the vectors  and  are combinations of those for electromagnetic and elastodynamic waves (rows <ref> and <ref>).

For seismoelectric waves (row <ref>), constitutive equations (<ref>) and (<ref>) account for the coupling.
In this case the vectors  and  are combinations of those for electromagnetic and poroelastodynamic waves (rows <ref> and <ref>).

In all cases, except for quantum mechanical waves, the vectors  and  are defined such that they consitute the power-flux density  in the x_3-direction via


    =(^†+^†),


where the superscript † denotes transposition and complex conjugation.
For quantum mechanical waves,  represents the probability current density in the x_3-direction.
Vectors _1 and _2 and operator matrices _11, _12, _21 and _22 
are defined in Appendices <ref> to <ref> for the different wave phenomena.
The operator matrices contain specific combinations of space-dependent medium parameters (or, for quantum mechanics, the potential V and mass m) 
and spatial differential operators ∂_1 and ∂_2
(standing for ∂/∂ x_1 and ∂/∂ x_2, respectively).

Equation (<ref>), with the operator matrix specified in the appendices, 
may be used as a starting point for generalising many of the applications mentioned in the introduction
(analysis of surface waves, derivation of long-wave effective medium parameters,
derivation of propagator matrices, decomposition into downgoing and upgoing waves, modelling wide-angle propagation in laterally variant media, 

etc.). A discussion of these applications is beyond the scope of this paper. 

Here we focus on the symmetry of the operator matrix and its use in unified reciprocity theorems and representations.




 §.§ Symmetry properties of the operator matrix


We discuss the symmetry properties of the operator matrix.
First, consider a scalar operator U, containing space-dependent parameters  and differential operators ∂_1 and ∂_2.
We introduce its transpose U^t and adjoint U^† via their integral properties


    ∫_ ( Uf)g  d^2=∫_ f( U^tg)  d^2


and


    ∫_ ( Uf)^*g  d^2=∫_ f^*( U^† g)  d^2.


Here  is the horizontal coordinate vector (x_1,x_2), superscript * denotes complex conjugation,
 denotes an infinite horizontal integration boundary at arbitrary depth x_3, and
f=f() and g=g() are space-dependent functions with sufficient decay along  towards infinity.
Equation (<ref>) implies 


    ( U V W)^t= W^t V^t U^t,


where also V and W are scalar operators.
Equations (<ref>) and (<ref>)  imply  


    U^†=( U^t)^*.


For the special case that U=∂_1, equation  (<ref>) implies
(via integration by parts) ∂_1^t=-∂_1.  Similarly,  ∂_2^t=-∂_2. Hence,


    ∂_α^t=-∂_α,


where Greek subscripts take on the values 1 and 2. 
Using this property and equation (<ref>), we find for example for the operator in equation (<ref>), 
(∂_α b_αβ∂_β)^t=∂_β b_αβ∂_α
(Einstein's summation convention applies to repeated subscripts). Since b_αβ=b_βα, this implies
(∂_α b_αβ∂_β)^t=∂_β b_βα∂_α=∂_α b_αβ∂_β and, using equation (<ref>), 
(∂_α b_αβ∂_β)^†=(∂_α b_αβ∂_β)^*=∂_α b_αβ^*∂_β. 



Next, we consider an operator matrix , of which the entries are operators containing space-dependent parameters and  differential operators ∂_1 and ∂_2.
Analogous to equations (<ref>) and (<ref>), we introduce its transpose ^t and its adjoint ^† via


    ∫_ ( f)^t g  d^2=∫_ f^t(^t g)  d^2


and


    ∫_ ( f)^† g  d^2=∫_ f^†(^† g)  d^2,


where f= f() and g= g() are space-dependent vector functions with sufficient decay  along  towards infinity.
Equation (<ref>) implies that  ^t involves transposition of the matrix and transposition of the operators contained in the matrix. 
For example, for a 2× 2 operator matrix , we have


    [ U_11 U_12; U_21 U_22 ]^t=[ U_11^t U_21^t; U_12^t U_22^t ].


Equation (<ref>) implies


    ()^t=^t^t^t,


where also  and  are operator matrices.
Equations (<ref>) and (<ref>) imply 


    ^†=(^t)^*.


Using equations (<ref>), (<ref>) and (<ref>), it follows that operator matrices _11, _12, _21 and _22, 
defined in Appendices <ref> to <ref>, obey the following symmetry relations


    _11^t    =    -_22,
    _12^t    =   _12, 
    _21^t    =   _21.



In Appendices  <ref> to <ref> we define adjoint medium parameters (or, for quantum mechanics, an adjoint potential).
When a medium is dissipative, its adjoint is effectual, and vice versa <cit.>.
A wave propagating through an effectual medium gains energy. 
Effectual media play a role in the reciprocity theorems and representations, discussed in sections <ref> and <ref>. 
An adjoint medium parameter is denoted by an overbar. An operator with an overbar means that the medium parameters  contained in that operator are replaced by their adjoints.
For example, for the operator in equation (<ref>) we have 
1/ıω∂_α b_αβ∂_β=1/ıω∂_αb̅_αβ∂_β.
Since b̅_αβ=b_αβ^* this becomes 
1/ıω∂_α b_αβ∂_β=1/ıω∂_α b_αβ^*∂_β=-(1/ıω∂_α b_αβ∂_β)^*.

For the operator matrices in Appendices  <ref> to <ref> in an adjoint medium we have


    _11   =   _11^*,
    _12   =   -_12^*,
    _21   =   -_21^*,
    _22   =   _22^*.


Using equation (<ref>), we find from equations (<ref>)  - (<ref>)


    _11^†   =    -_22,
    _12^†   =   -_12,
    _21^†   =    -_21,
    _22^†   =    -_11.


From equations (<ref>)  - (<ref>), we find for the 
operator matrix   defined in equation (<ref>)


    ^t   =   -,
    ^*   =   ,
    ^†   =   -,


with


    =[   O   I; - I   O ],
      =[   I   O;   O - I ],
      =[ O I; I O ],


where O and I are zero and identity matrices of appropriate size.

Symmetry relations in the wavenumber-frequency domain for the special case of a laterally invariant medium (or potential) are given in Appendix <ref>.



§ MATRIX-VECTOR WAVE FIELD RECIPROCITY THEOREMS


In wave theory, a reciprocity theorem formulates a mathematical relation between two 
states (wave fields, sources and medium parameters) 
in the same spatial domain.
An early reference for the acoustic reciprocity theorem is <cit.>, who referred to it as Helmholtz's theorem.
<cit.> formulated a reciprocity theorem for electromagnetic fields.
Early references for elastodynamic reciprocity theorems are <cit.> and <cit.>. 
<cit.> and <cit.> formulated reciprocity theorems for piezoelectric and seismoelectric waves, respectively. 
Comprehensive overviews of the history
of reciprocity theorems  and their applications are given by <cit.>, <cit.> and <cit.>.



Matrix-vector wave equation (<ref>) and symmetry relations (<ref>) and (<ref>) underly unified matrix-vector reciprocity theorems 
in an inhomogeneous medium (or potential). We consider two 
states A and B, characterised by independent wave  vectors _A(,ω) and _B(,ω), obeying matrix-vector wave equation (<ref>),
with source vectors _A(,ω) and _B(,ω), and operator matrices _A(,ω) and _B(,ω).
The subscripts A and B of these operator matrices refer to the, possibly different, medium parameters in states A and B.
We assume that outside a finite domain, the medium (or potential) and its adjoint are lossless in both states.
We consider a spatial domain  enclosed by two infinite horizontal boundaries _0 and _1 (with _1 below _0), together denoted by , see Figure <ref>.
The starting point for deriving reciprocity theorems for the wave fields in states A and B is given by
the quantities ∂_3{_A^t_B} and  ∂_3{_A^†_B} in domain . 
Applying the product rule for differentiation gives


    ∂_3{_A^t_B}   =    (∂_3_A^t)_B + _A^t(∂_3_B),
    ∂_3{_A^†_B}   =    (∂_3_A^†)_B + _A^†(∂_3_B).


Note that equation (<ref>), which defines the transpose of the horizontal differential operator ∂_α, 

does not apply to the vertical differential operator ∂_3. 
Hence, we may replace (∂_3_A^t) by (∂_3_A)^t in equation (<ref>), and (∂_3_A^†) by (∂_3_A)^† in equation (<ref>).
Using wave equation (<ref>) for both states in the right-hand sides of equations (<ref>)  and (<ref>),
integrating both sides of these equations over domain  and applying the theorem of Gauss to the left-hand sides, we obtain 


    ∫__A^t_B n_3 d^2=∫_[((_A_A)^t+_A^t)_B +_A^t(_B_B+_B) ] d^3


and


    ∫__A^†_B n_3 d^2=∫_[((_A_A)^†+_A^†)_B +_A^†(_B_B+_B) ] d^3.


Here  n_3 is the vertical component of the outward pointing normal vector on , with n_3=-1 at the upper boundary _0
and n_3=+1 at the lower boundary _1. 
The integrals on the right-hand sides can be written as 


    ∫_{⋯} d^3=∫_x_3,0^x_3,1 dx_3∫_{⋯} d^2,


where x_3,0 and x_3,1 denote the depths of _0 and _1, respectively.
Hence, at each depth level between _0 and _1 we can use the integral properties of transpose and adjoint operators, as formulated by
 equations (<ref>) and (<ref>). Together with the 
symmetry relations  (<ref>) and (<ref>) for operator _A, we thus obtain the following matrix-vector reciprocity theorems 


    ∫_(_A^t_B +_A^t_B ) d^3=
    ∫__A^t_B n_3 d^2+∫__A^t(_A-_B)_B  d^3


and


    ∫_(_A^†_B + _A^†_B) d^3=
    ∫__A^†_B n_3 d^2+∫__A^†(_A-_B)_B  d^3.


Equation (<ref>) is a convolution-type reciprocity theorem <cit.> because products like _A^t_B
in the frequency domain correspond to convolutions in the time domain. Equation (<ref>) is a correlation-type reciprocity theorem <cit.> 
because products like _A^†_B in the frequency domain correspond to correlations in the time domain. 
These matrix-vector reciprocity  theorems have been previously derived for acoustic waves <cit.>.
Because these theorems follow from the unified matrix-vector wave equation (<ref>), with
unified symmetry relations (<ref>) and (<ref>), they hold for all wave phenomena listed in Table <ref>.
In the next section we use these theorems as the basis for matrix-vector wave field representations. Here we consider some special cases
of these theorems.


- Power balance

When the sources, medium parameters and wave fields are identical in both states, we may drop the subscripts A and B. In this case equation (<ref>) simplifies to


    ∫_(^† + ^†) d^3=
    ∫_^† n_3 d^2+∫_^†(-) d^3.


Because ^†=(^†+^†)=,
the first term on the right-hand side is the power flux (or probability current) through the  boundary =_0∪_1 (i.e., the power leaving the domain ). 
Hence, equation (<ref>) formulates the unified power  balance.
The term on the left-hand side is the power generated by the sources in  and
the second term on the right-hand side is the dissipated power in .


- Propagation invariants

When there are no sources in  and the medium parameters in  are equal in the two states, the domain integrals in equation (<ref>) vanish, hence



    ∫__0∪_1_A^t_B n_3 d^2=0,


or, since n_3=-1 at  _0 and n_3=+1 at  _1,


    ∫__0_A^t_B d^2 = ∫__1_A^t_B d^2.


Since this holds for any choice of the domain , we infer that the quantity 


    ∫__A^t_B  d^2,


with  denoting a horizontal plane at arbitrary depth x_3,
is a unified propagation invariant (i.e., it is independent of the depth x_3 of ).
On the other hand, when the medium parameters 

are each other's adjoints in the two states, 
we find in a similar way from equation (<ref>) that the quantity



    ∫__A^†_B  d^2


is a unified propagation invariant. 
Propagation invariants have been extensively used in the 
analysis of symmetry properties of reflection and transmission responses and for the design of efficient 
numerical modelling schemes for acoustic and elastodynamic wave fields <cit.>.



§ MATRIX-VECTOR WAVE FIELD REPRESENTATIONS




 §.§ Representation of the convolution type


A wave field representation is obtained by replacing one of the states in a reciprocity theorem by a Green's state 

<cit.>. Here we derive a unified matrix-vector wave field representation  

from the matrix-vector reciprocity theorem of  the convolution type (equation <ref>).


We introduce the Green's matrix (,_A,ω) 
(with the same dimensions as matrix )
as the solution of the unified matrix-vector wave equation (<ref>), with the source vector 
 replaced by a diagonal point-source matrix. Hence


    ∂_3 = + Iδ(-_A),


where I is an identity matrix and _A defines the position of the point source.
We let  represent the forward propagating solution of equation (<ref>), 
which corresponds to imposing causality in the time domain, i.e., (,_A,t)= O for t<0, where O is a zero matrix
(the relation between functions in the time- and frequency domain is defined by  the Fourier transform, equation (<ref>)).

Before we derive the unified wave field representation, we first derive a reciprocity relation for the Green's matrix. To this end
we define a second forward propagating Green's matrix (,_B,ω), with its point source at _B.
We assume that _A and _B are both situated in . 

We replace   _A and _B in reciprocity theorem (<ref>) by (,_A,ω) and
(,_B,ω), respectively. Accordingly, we replace  _A and _B by Iδ(-_A) and Iδ(-_B), respectively.

Both Green's matrices are defined in the same medium, hence, _A=_B. This implies that 
the second integral on the right-hand side of equation (<ref>) vanishes.

When Neumann or Dirichlet boundary conditions apply on , 
or when the medium outside  is homogeneous, the first integral on the right-hand side of equation (<ref>) vanishes as well.
We thus obtain


    (_A,_B,ω) +^t(_B,_A,ω)= O.


Using ^-1=- this gives



    (_A,_B,ω)=^t(_B,_A,ω),


which is the unified source-receiver reciprocity relation for the Green's matrix. 
 
Next, we use the same reciprocity theorem to derive a representation for the actual wave field vector . 

We let state B be the 
actual state (i.e., actual wave field, source and medium parameters). For convenience we drop the subscript B from _B, _B and _B. 
For state A we choose again the Green's state.


Hence, 
 we replace _A by (,_A,ω) and  
_A by Iδ(-_A).
Operator _A may be defined in a reference medium or in the actual medium.


Making these substitutions in equation (<ref>), pre-multiplying all terms by -, using  equation  (<ref>) and -= I, we obtain


    (_A)(_A,ω)   =   ∫_(_A,,ω)(,ω) d^3
    -∫_(_A,,ω)(,ω)n_3 d^2
       +   ∫_(_A,,ω){-_A}(,ω) d^3,


where (_A) is the characteristic function, defined as


    (_A)=
    
    1,      for _A inside , 
    ,    for _A on ,
    
    0,     for _A outside .


The left-hand side of equation (<ref>)  is the actual wave field vector , observed at _A (assuming _A is inside ). 
The right-hand side contains, respectively, a contribution from the source distribution (,ω) inside , 
a contribution from the wave field (,ω) at the boundary , and
a contribution from the contrast operator -_A, applied to the wave field (,ω) inside . 
This unified matrix-vector wave field representation holds for all wave phenomena listed in Table <ref>.

This representation can often be simplified, which leads to different applications.
For example, when the medium outside the domain   is homogeneous, source free and identical in both states, 
the boundary integral on the right-hand side vanishes. 
The remaining representation (with _A defined in a reference medium) forms a basis for the analysis of forward scattering problems.
On the other hand, when _A is defined in the actual medium (i.e., _A=) and the domain   is source free, only the boundary integral on the right-hand side remains.
In this case, equation (<ref>)
is a generalisation of the Kirchhoff-Helmholtz integral <cit.>, which finds applications 
in forward wave field extrapolation problems.



 §.§ Representation of the correlation type


Representations of the correlation type find their application in inverse source problems <cit.>, inverse scattering problems
<cit.>, 
imaging <cit.>, time-reversal acoustics
<cit.>, and Green's function retrieval from ambient noise <cit.>.
There are several ways to approach the representation of the correlation type.
The homogeneous Green's function representation <cit.> elegantly covers most of the aforementioned applications for 
scalar wave fields. It is obtained by replacing both states in the reciprocity theorem of the correlation type by Green's states.
Here we derive a unified representation for the homogeneous Green's matrix by substituting
two Green's matrices into the matrix-vector reciprocity theorem of the correlation type (equation <ref>).

Before we discuss the homogeneous Green's matrix, we introduce the Green's matrix of the adjoint medium, (,_A,ω), 
as the forward propagating solution of the  following matrix-vector wave equation


    ∂_3=+ Iδ(-_A).


Pre- and post multiplying all terms by  and subsequently using equation (<ref>) gives 


    ∂_3 =^* +δ(-_A).


Taking the complex conjugate of all terms and using = I gives


    ∂_3^* =^* + Iδ(-_A).


Subtracting all terms in this equation from the corresponding terms in equation (<ref>) we obtain 


    ∂_3_ h(,_A,ω) =_ h(,_A,ω),


with


    _ h(,_A,ω)=(,_A,ω)-^*(,_A,ω).


Because _ h(,_A,ω) obeys a matrix-vector wave equation without a source term, we call it the homogeneous Green's matrix.
The second term on the right-hand side represents a backward propagating wave field in the adjoint medium.

Next, we use the correlation-type reciprocity theorem (equation (<ref>)) to derive a representation for the  homogeneous Green's matrix _ h.
For state A we choose the Green's matrix in the adjoint medium medium, hence, we replace _A by (,_A,ω), _A by Iδ(-_A), and _A by .
For state B we choose the Green's matrix in the actual medium, hence, we replace _B by (,_B,ω), _B by Iδ(-_B), and _B by .
With these choices the contrast operator _A-_B=    - vanishes.
Making these substitutions in equation (<ref>), taking _A and _B both inside ,
pre-multiplying all terms by , using  equations  (<ref>) and (<ref>), = I and ==-, we obtain


    _ h(_A,_B,ω)=∫_^†(,_A,ω)(,_B,ω)n_3 d^2.


This unified homogeneous Green's matrix representation holds for all wave phenomena listed in Table <ref>.
It forms the basis for  generalising the applications mentioned at the beginning of this section.



§ CONCLUSIONS


A unified  matrix-vector wave equation is presented for acoustic, quantum mechanical, electromagnetic, elastodynamic, 
poroelastodynamic, piezoelectric and seismoelectric waves. For most cases a 3D inhomogeneous, anisotropic, dissipative medium is considered. 
The unified equation may be used as a basis for generalising various applications of the elastodynamic matrix-vector wave equation, such as
the analysis of surface waves, the derivation of long-wave effective medium parameters,
the derivation of propagator matrices, decomposition into downgoing and upgoing waves, modelling wide-angle propagation in laterally variant media, etc.

The operator matrix in the matrix-vector wave equation obeys unified  symmetry relations.
These symmetry relations underly unified reciprocity theorems of the convolution and correlation type, which, in turn, 
form the basis for representations of the wave vector and the homogeneous Green's matrix.
Reciprocity theorems and representations find applications in
forward modelling problems, inverse source and inverse scattering problems, imaging, time-reversal methods and Green's function retrieval from ambient noise.
The unified treatment in this paper provides a starting point for generalising these applications to a broad range of wave phenomena.





§ ACKNOWLEDGMENTS


The author thanks Evert Slob and Niels Grobbe for fruitful discussions. The constructive comments of the  reviewers and editors
(including Richard Gibson, Rune Mittet and Jörg Renner)
are very much appreciated and helped a lot to improve the paper. 
This research  has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research 
and innovation programme (grant agreement No: 742703).

gji


93
natexlab#1#1

[Achenbach(1973)]Achenbach73Book
Achenbach, J. D., 1973.
Wave propagation in elastic solids, North Holland.

[Achenbach(2003)]Achenbach2003Book
Achenbach, J. D., 2003.
Reciprocity in elastodynamics, Cambridge University Press.

[Aki & Richards(1980)]Aki80Book
Aki, K. & Richards, P. G., 1980.
Quantitative seismology, Vol. I, W.H. Freeman and
  Company, San Fransisco.

[Auld(1973)]Auld73Book
Auld, B. A., 1973.
Acoustic fields and waves in solids, Wiley-Interscience, New
  York.

[Auld(1979)]Auld79WM
Auld, B. A., 1979.
General electromechanical reciprocity relations applied to the
  calculation of elastic wave scattering coefficients, Wave Motion,
  1, 3–10.

[Backus(1962)]Backus62JGR
Backus, G. E., 1962.
Long-wave elastic anisotropy produced by horizontal layering, J. Geophys. Res., 67(11), 4427–4440.

[Berkhout(1982)]Berkhout82Book
Berkhout, A. J., 1982.
Seismic Migration. Imaging of acoustic energy by wave field
  extrapolation. A. Theoretical aspects, Elsevier, Amsterdam.

[Berkhout(1985)]Berkhout85Book
Berkhout, A. J., 1985.
Seismic Migration. Imaging of acoustic energy by wave field
  extrapolation. A. Theoretical aspects (3rd edition), Elsevier,
  Amsterdam.

[Biot(1956a)]Biot56JASA
Biot, M. A., 1956a.
Theory of propagation of elastic waves in a fluid-saturated porous
  solid: I. Low frequency range, J. Acoust. Soc. Am., 28, 168–178.

[Biot(1956b)]Biot56JASA2
Biot, M. A., 1956b.
Theory of propagation of elastic waves in a fluid-saturated porous
  solid: II. High frequency range, J. Acoust. Soc. Am., 28, 179–191.

[Bleistein(1984)]Bleistein84Book
Bleistein, N., 1984.
Mathematical methods for wave phenomena, Academic Press,
  Inc., Orlando.

[Bojarski(1983)]Bojarski83JASA
Bojarski, N. N., 1983.
Generalized reaction principles and reciprocity theorems for the wave
  equations, and the relationship between the time-advanced and time-retarded
  fields, J. Acoust. Soc. Am., 74,
  281–285.

[Born & Wolf(1965)]Born65Book
Born, M. & Wolf, E., 1965.
Principles of optics, Pergamon Press, London.

[Carcione(2007)]Carcione2007Book
Carcione, J. M., 2007.
Wave fields in real media, Elsevier, Amsterdam.

[Chapman(1994)]Chapman94GJI
Chapman, C. H., 1994.
Reflection/transmission coefficients reciprocities in anisotropic
  media, Geophys. J. Int., 116, 498–501.

[Corones(1975)]Corones75JMAA
Corones, J. P., 1975.
Bremmer series that correct parabolic approximations, J. Math.  Anal. Appl., 50, 361–372.

[Dahlen & Tromp(1998)]Dahlen98Book
Dahlen, F. A. & Tromp, J., 1998.
Theoretical global seismology, Princeton University
  Press.

[de Hoop(1958)]Hoop58PHD
de Hoop, A. T., 1958.
Representation theorems for the displacement in an elastic solid
  and their applications to elastodynamic diffraction theory, Ph.D. thesis,
  Delft University of Technology, Delft (repository.tudelft.nl).

[de Hoop(1966)]Hoop66ASR
de Hoop, A. T., 1966.
An elastodynamic reciprocity theorem for linear, viscoelastic media,
  Appl. Sci. Res., 16, 39–45.

[de Hoop(1987)]Hoop87RS
de Hoop, A. T., 1987.
Time-domain reciprocity theorems for electromagnetic fields in
  dispersive media, Radio Sci., 22(7), 1171–1178.

[de Hoop(1988)]Hoop88JASA
de Hoop, A. T., 1988.
Time-domain reciprocity theorems for acoustic wave fields in fluids
  with relaxation, J. Acoust. Soc. Am., 84, 1877–1882.

[de Hoop(1995)]Hoop95Book
de Hoop, A. T., 1995.
Handbook of radiation and scattering of waves, Academic
  Press, London.

[de Hoop(1996a)]Hoop96GJI
de Hoop, A. T., 1996a.
A general correspondence principle for time-domain electromagnetic
  wave and diffusion fields, Geophys. J. Int., 127, 757–761.

[de Hoop(1996b)]Hoop96JMP
de Hoop, M. V., 1996b.
Generalization of the Bremmer coupling series, J. Math.  Phys., 37, 3246–3282.

[Derode et al.(2003)Derode, Larose, Tanter, de Rosny, Tourin,
  Campillo, & Fink]Derode2003JASA
Derode, A., Larose, E., Tanter, M., de Rosny, J., Tourin, A., Campillo, M.,
  & Fink, M., 2003.
Recovering the Green's function from field-field correlations in an
  open scattering medium (L), J. Acoust. Soc. Am., 113(6), 2973–2976.

[Devaney(1982)]Devaney82UI
Devaney, A. J., 1982.
A filtered backpropagation algorithm for diffraction tomography, Ultrasonic Imaging, 4, 336–350.

[Esmersoy & Oristaglio(1988)]Esmersoy88GEO
Esmersoy, C. & Oristaglio, M., 1988.
Reverse-time wave-field extrapolation, imaging, and inversion, Geophysics, 53, 920–931.

[Fink & Prada(2001)]Fink2001IP
Fink, M. & Prada, C., 2001.
Acoustic time-reversal mirrors, Inverse Probl., 17,
  R1–R38.

[Fishman(1992)]Fishman92JMP
Fishman, L., 1992.
Exact and operator rational approximate solutions of the Helmholtz,
  Weyl composition equation in underwater acoustics - The quadratic
  profile, J. Math. Phys., 33(5), 1887–1914.

[Fishman & McCoy(1984)]Fishman84JMP
Fishman, L. & McCoy, J. J., 1984.
Derivation and application of extended parabolic wave theories. I.
  The factorized Helmholtz equation, J. Math. Phys., 25(2),
  285–296.

[Fokkema & van den Berg(1993)]Fokkema93Book
Fokkema, J. T. & van den Berg, P. M., 1993.
Seismic applications of acoustic reciprocity, Elsevier,
  Amsterdam.

[Frasier(1970)]Frasier70GEO
Frasier, C. W., 1970.
Discrete time solution of plane P-SV waves in a plane layered
  medium, Geophysics, 35, 197–219.

[Frazer & Fryer(1989)]Frazer89GJI
Frazer, L. N. & Fryer, G. J., 1989.
Useful properties of the system matrix for a homogeneous anisotropic
  visco-elastic solid, Geophys. J. Int., 97,
  173–177.

[Frazer & Sen(1985)]Frazer85GJRAS
Frazer, L. N. & Sen, M. K., 1985.
Kirchhoff-Helmholtz reflection seismograms in a laterally
  inhomogeneous multi-layered elastic medium, I, Theory, Geophys. J. R. Astr. Soc., 80, 121–147.

[Gangi(1970)]Gangi70JGR
Gangi, A. F., 1970.
A derivation of the seismic representation theorem using seismic
  reciprocity, J. Geophys. Res., 75, 2088–2095.

[Gangi(2000)]Gangi2000GJI
Gangi, A. F., 2000.
Constitutive equations and reciprocity, Geophys. J. Int., 143, 311–318.

[Gelinsky & Shapiro(1997)]Gelinsky97GEO
Gelinsky, S. & Shapiro, S. A., 1997.
Poroelastic Backus averaging for anisotropic layered fluid- and
  gas-saturated sediments, Geophysics, 62(6), 1867–1878.

[Gilbert & Backus(1966)]Gilbert66GEO
Gilbert, F. & Backus, G. E., 1966.
Propagator matrices in elastic wave and vibration problems, Geophysics, 31(2), 326–332.

[Grobbe(2016)]Grobbe2016PHD
Grobbe, N., 2016.
Coupled poroelastic waves and electromagnetic fields in layered
  media: Theory, modeling, and interferometric synthesis, Ph.D. thesis,
  Delft University of Technology, Delft (repository.tudelft.nl).

[Gurevich & Lopatnikov(1995)]Gurevich95GJI
Gurevich, B. & Lopatnikov, S. L., 1995.
Velocity and attenuation of elastic waves in finely layered porous
  rocks, Geophys. J. Int., 121, 933–947.

[Haartsen & Pride(1997)]Haartsen97JGR
Haartsen, M. W. & Pride, S. R., 1997.
Electroseismic waves from point sources in layered media, J. Geophys. Res., 102, 24745–24769.

[Haines(1988)]Haines88GJI
Haines, A. J., 1988.
Multi-source, multi-receiver synthetic seismograms for laterally
  heterogeneous media using F-K domain propagators, Geophys. J. Int., 95, 237–260.

[Haines & de Hoop(1996)]Haines96JMP
Haines, A. J. & de Hoop, M. V., 1996.
An invariant imbedding analysis of general wave scattering problems,
  J. Math. Phys., 37, 3854–3881.

[Haskell(1953)]Haskell53BSSA
Haskell, N. A., 1953.
The dispersion of surface waves on multilayered media, Bull. Seism. Soc. Am., 43, 17–34.

[Honein et al.(1991)Honein, Braga, Barbone, & Herrmann]Honein91JIMSS
Honein, B., Braga, A. M. B., Barbone, P., & Herrmann, G., 1991.
Wave propagation in piezoelectric layered media with some
  applications, J. Intel. Mat. Syst. Struct.,
  2, 542–557.

[Kennett(1983)]Kennett83Book
Kennett, B. L. N., 1983.
Seismic wave propagation in stratified media, Cambridge
  University Press.

[Kennett & Illingworth(1981)]Kennett81GJRAS
Kennett, B. L. N. & Illingworth, M. R., 1981.
Seismic waves in a stratified half-space - III, Piecewise smooth
  models, Geophys. J. R. Astr. Soc., 66, 633–675.

[Kennett & Kerry(1979)]Kennett79GJRAS
Kennett, B. L. N. & Kerry, N. J., 1979.
Seismic waves in a stratified half-space, Geophys. J. R. Astr. Soc., 57, 557–584.

[Kennett et al.(1978)Kennett, Kerry, & Woodhouse]Kennett78GJRAS
Kennett, B. L. N., Kerry, N. J., & Woodhouse, J. H., 1978.
Symmetries in the reflection and transmission of elastic waves, Geophys. J. R. Astr. Soc., 52, 215–230.

[Kennett et al.(1990)Kennett, Koketsu, & Haines]Kennett90GJI
Kennett, B. L. N., Koketsu, K., & Haines, A. J., 1990.
Propagation invariants, reflection and transmission in anisotropic,
  laterally heterogeneous media, Geophys. J. Int.,
  103, 95–101.

[Knopoff(1956)]Knopoff56JASA
Knopoff, L., 1956.
Diffraction of elastic waves, J. Acoust. Soc. Am., 28, 217–229.

[Knopoff & Gangi(1959)]Knopoff59GEO
Knopoff, L. & Gangi, A. F., 1959.
Seismic reciprocity, Geophysics, 24, 681–691.

[Koketsu et al.(1991)Koketsu, Kennett, & Takenaka]Koketsu91GJI
Koketsu, K., Kennett, B. L. N., & Takenaka, H., 1991.
2-D reflectivity method and synthetic seismograms for irregularly
  layered structures - II. Invariant embedding approach, Geophys. J. Int., 105, 119–130.

[Landau & Lifshitz(1960)]Landau60Book
Landau, L. D. & Lifshitz, E. M., 1960.
Electrodynamics of continuous media, Pergamon Press.

[Lindsey & Braun(2004)]Lindsey2004AJSS
Lindsey, C. & Braun, D. C., 2004.
Principles of seismic holography for diagnostics of the shallow
  subphotosphere, Astrophys. J. Suppl. Series, 155(1), 209–225.

[Lorentz(1895)]Lorentz1895KNAW
Lorentz, H. A., 1895.
The theorem of poynting concerning the energy in the electromagnetic
  field and two general propositions concerning the propagation of light, Verslagen der Afdeeling Natuurkunde van de Koninklijke Akademie van
  Wetenschappen, 4, 176–187.

[Løseth & Ursin(2007)]Loseth2007GJI
Løseth, L. O. & Ursin, B., 2007.
Electromagnetic fields in planarly layered anisotropic media, Geophys. J. Int., 170, 44–80.

[Mavko et al.(2009)Mavko, Mukerji, & Dvorkin]Mavko2009Book
Mavko, G., Mukerji, T., & Dvorkin, J., 2009.
The rock physics handbook: Tools for seismic analysis of
  porous media, Cambridge University Press.

[Maynard et al.(1985)Maynard, Williams, & Lee]Maynard85JASA
Maynard, J. D., Williams, E. G., & Lee, Y., 1985.
Nearfield acoustic holography: I. Theory of generalized
  holography and the development of NAH, J. Acoust. Soc. Am., 78(4), 1395–1413.

[Merzbacher(1961)]Merzbacher61Book
Merzbacher, E., 1961.
Quantum mechanics, John Wiley and Sons, Inc., New York.

[Messiah(1961)]Messiah61Book
Messiah, A., 1961.
Quantum mechanics, Volume I, North-Holland Publishing
  Company, Amsterdam.

[Mittet(2015)]Mittet2015GEO
Mittet, R., 2015.
Seismic wave propagation concepts applied to the interpretation of
  marine controlled-source electromagnetics, Geophysics, 80(2),
  E63–E81.

[Morse & Feshbach(1953)]Morse53Book
Morse, P. M. & Feshbach, H., 1953.
Methods of theoretical physics, Vol. I, McGraw-Hill
  Book Company Inc., New York.

[Müller(1985)]Muller85JG
Müller, G., 1985.
The reflectivity method: a tutorial, J. Geophys., 58(1-3), 153–174.

[Norris(1993)]Norris93JASA
Norris, A. N., 1993.
Low-frequency dispersion and attenuation in partially saturated
  rocks, J. Acoust. Soc. Am., 94(1),
  359–370.

[Oristaglio(1989)]Oristaglio89IP
Oristaglio, M. L., 1989.
An inverse scattering formula that uses all the data, Inverse Probl., 5, 1097–1105.

[Pao & Varatharajulu(1976)]Pao76JASA
Pao, Y. H. & Varatharajulu, V., 1976.
Huygens' principle, radiation conditions, and integral formulations
  for the scattering of elastic waves, J. Acoust. Soc. Am., 59, 1361–1371.

[Porter(1970)]Porter70JOSA
Porter, R. P., 1970.
Diffraction-limited, scalar image formation with holograms of
  arbitrary shape, J. Opt. Soc.  Am., 60,
  1051–1059.

[Porter & Devaney(1982)]Porter82JOSA
Porter, R. P. & Devaney, A. J., 1982.
Holography and the inverse source problem, J. Opt. Soc.  Am., 72, 327–330.

[Pride(1994)]Pride94PRB
Pride, S., 1994.
Governing equations for the coupled electromagnetics and acoustics of
  porous media, Phys. Rev. B, 50, 15678–15696.

[Pride & Haartsen(1996)]Pride96JASA
Pride, S. R. & Haartsen, M. W., 1996.
Electroseismic wave properties, J. Acoust. Soc. Am., 100, 1301–1315.

[Pride et al.(1992)Pride, Gangi, & Morgan]Pride92JASA
Pride, S. R., Gangi, A. F., & Morgan, F. D., 1992.
Deriving the equations of motion for porous isotropic media, J. Acoust. Soc. Am., 92(6), 3278 – 3290.

[Rayleigh(1878)]Rayleigh78Book
Rayleigh, J. W. S., 1878.
The theory of sound. Volume II, Dover Publications,
  Inc. (Reprint 1945).

[Schneider(1978)]Schneider78GEO
Schneider, W. A., 1978.
Integral formulation for migration in two and three dimensions, Geophysics, 43, 49–76.

[Schoenberg & Sen(1983)]Schoenberg83JASA
Schoenberg, M. & Sen, P. N., 1983.
Properties of a periodically stratified acoustic half-space and its
  relation to a Biot fluid, J. Acoust. Soc. Am., 73, 61–67.

[Takenaka et al.(1993)Takenaka, Kennett, & Koketsu]Takenaka93WM
Takenaka, H., Kennett, B. L. N., & Koketsu, K., 1993.
The integral operator representation of propagation invariants for
  elastic waves in irregularly layered media, Wave Motion, 17,
  299–317.

[Thomson(2015a)]Thomson2015GJI1
Thomson, C. J., 2015a.
Seismic transmission operator reciprocity - I: a direct derivation
  showing the role of elastic impedance operator symmetry, Geophys. J. Int., 202, 1064–1070.

[Thomson(2015b)]Thomson2015GJI2
Thomson, C. J., 2015b.
Seismic transmission operator reciprocity - II: impedance-operator
  symmetry via elastic lateral modes, Geophys. J. Int., 202, 1071–1087.

[Thomson(1950)]Thomson50JAP
Thomson, W. T., 1950.
Transmission of elastic waves through a stratified solid medium, J. Appl. Phys., 21, 89–93.

[Ursin(1983)]Ursin83GEO
Ursin, B., 1983.
Review of elastic and electromagnetic wave propagation in
  horizontally layered media, Geophysics, 48, 1063–1081.

[van Stralen(1997)]Stralen97PHD
van Stralen, M. J. N., 1997.
Directional decomposition of electromagnetic and acoustic
  wave-fields - Applications in integrated optics, exploration seismics and underwater acoustics, Ph.D. thesis, Delft University of Technology, Delft (repository.tudelft.nl).

[Wang & Rokhlin(2002)]Wang2002APL
Wang, L. & Rokhlin, S. I., 2002.
Recursive asymptotic stiffness matrix method for analysis of surface
  acoustic wave devices on layered piezoelectric media, Appl. Phys.
  Lett., 81(21), 4049–4051.

[Wapenaar(1996a)]Wapenaar96GJI2
Wapenaar, C. P. A., 1996a.
One-way representations of seismic data, Geophys. J. Int., 127, 178–188.

[Wapenaar(1996b)]Wapenaar96JASA
Wapenaar, C. P. A., 1996b.
Reciprocity theorems for two-way and one-way wave vectors: a
  comparison, J. Acoust. Soc. Am., 100,
  3508–3518.

[Wapenaar & Berkhout(1989)]Wapenaar89Book
Wapenaar, C. P. A. & Berkhout, A. J., 1989.
Elastic wave field extrapolation, Elsevier, Amsterdam.

[Wapenaar & Grimbergen(1996)]Wapenaar96GJI1
Wapenaar, C. P. A. & Grimbergen, J. L. T., 1996.
Reciprocity theorems for one-way wave fields, Geophys. J. Int., 127, 169–177.

[Wapenaar et al.(2001)Wapenaar, Dillen, & Fokkema]Wapenaar2001RS
Wapenaar, C. P. A., Dillen, M. W. P., & Fokkema, J. T., 2001.
Reciprocity theorems for electromagnetic or acoustic one-way wave
  fields in dissipative inhomogeneous media, Radio Sci., 36,
  851–863.

[Wapenaar(2003)]Wapenaar2003GEO
Wapenaar, K., 2003.
Synthesis of an inhomogeneous medium from its acoustic transmission
  response, Geophysics, 68, 1756–1759.

[Weaver & Lobkis(2004)]Weaver2004JASA
Weaver, R. L. & Lobkis, O. I., 2004.
Diffuse fields in open systems and the emergence of the Green's
  function (L), J. Acoust. Soc. Am., 116(5), 2731–2734.

[Weston(1989)]Weston89JMP
Weston, V. H., 1989.
Wave splitting and the reflection operator for the wave equation in
  ℝ^3, J. Math. Phys., 30, 2545–2562.

[White & Zhou(2006)]White2006SIAM
White, B. S. & Zhou, M., 2006.
Electroseismic prospecting in layered media, SIAM J. Appl.  Math., 67, 69–98.

[Willis(2012)]Willis2012CRM
Willis, J. R., 2012.
The construction of effective relations for waves in a composite,
  Comptes Rendus Mecanique, 340, 181–192.

[Woodhouse(1974)]Woodhouse74GJR
Woodhouse, J. H., 1974.
Surface waves in a laterally varying layered structure, Geophys. J. R. Astr. Soc., 37, 461–490.

[Zhao et al.(2012)Zhao, Pan, & Zhong]Zhao2012JAP
Zhao, J., Pan, Y., & Zhong, Z., 2012.
Theoretical study of shear horizontal wave propagation in
  periodically layered piezoelectric structure, J. Appl.
  Phys., 111, 064906.









§ ACOUSTIC WAVES


The basic equations for acoustic wave propagation 

are  the linearised  equation of motion


    ∂_t _i+∂_i =f_i


and the linearised deformation equation


    -∂_t Θ + ∂_i v_i=q


<cit.>. Here _i=_i(,t) is the momentum density as a function of spatial position  and time t,
 =(,t) is the acoustic pressure,
  Θ=Θ(,t) the cubic dilatation,
 v_i=v_i(,t) the particle velocity, and 
 f_i=f_i(,t) and q=q(,t)  represent the sources in terms of external force density and volume-injection rate density, respectively (the function q should not be confused with vector  and its components  and   in equations (<ref>) and (<ref>)).
Operator ∂_i stands for differentiation in the x_i-direction. 
Lower-case Latin subscripts (except t) take on the values 1, 2 and 3, and Einstein's summation convention applies to repeated subscripts.
Operator ∂_t  stands for the temporal differential operator ∂/∂ t.
The constitutive relations for an inhomogeneous, anisotropic fluid are given by


    _i   =   ρ_ijv_j,
    Θ   =   -κ,


where ρ_ij=ρ_ij() and κ=κ()  are the mass density and  compressibility, respectively.
To account for anisotropy, the mass density is defined as a tensor. Although ideal fluids are by definition isotropic, inhomogeneities at the micro 
scale can often be represented by effective anisotropic parameters at the macro scale.
For example, a periodic stratified fluid can, in the long wavelength limit, be represented by a homogeneous fluid with an effective transverse isotropic mass density tensor 
and an effective isotropic compressibility <cit.>. The mass density tensor is symmetric, i.e., ρ_ij=ρ_ji.
Substituting the constitutive relations (<ref>) and (<ref>) into equations (<ref>) and (<ref>) yields


    ρ_ij∂_t v_j+∂_i    =   f_i,
    κ∂_t  + ∂_i v_i   =   q.





We define the temporal Fourier transform of a space- and time-dependent function h(,t) as


    h(,ω)=∫_-∞^∞ h(,t) exp(ıω t) dt,


where ı is the imaginary unit.
For notational convenience, we use the same symbol (here h) for quantities in the time domain and in the frequency domain.
We use equation  (<ref>) to transform equations  (<ref>) and (<ref>)
to the frequency domain. The time derivatives are thus replaced by -ıω,
hence


    -ıωρ_ij v_j + ∂_i p   =   f_i,
    
    -ıωκ p + ∂_i v_i   =   q,


with  p=p(,ω), v_i=v_i(,ω),   f_i=f_i(,ω) and q=q(,ω). 
In a lossless medium, the parameters  ρ_ij() and κ() are real-valued and frequency independent.
To account for losses, we replace them by complex-valued, frequency-dependent parameters ρ_ij=ρ_ij(,ω) and
κ=κ(,ω) <cit.>.

The quantities p and v_3 constitute the power-flux density  in the x_3-direction, via


    =(p^*v_3+v_3^*p).


We choose these quantities for the 1× 1 vectors  and  in equations (<ref>) and (<ref>), hence


    =p,  =v_3.


To arrive at a set of equations for these quantities, 
we need to eliminate the remaining wave field quantities, v_1 and v_2, from equations (<ref>) and (<ref>).
To this end, we first introduce the inverse of the mass density tensor, the so-called specific volume tensor _hi=_hi(,ω), via


    _hiρ_ij=δ_hj,


with δ_hj denoting the Kronecker delta. On account of the symmetry of the mass density tensor and equation (<ref>), 
the specific volume tensor is symmetric as well,  hence _hi=_ih.
Applying _hi to both sides of equation (<ref>), using equation (<ref>),  gives


    -ıω v_h + _hi∂_i p=_hif_i.


We separate the derivatives in the x_3-direction from the lateral derivatives in equations (<ref>) and (<ref>), according to


    ∂_3 p   =   _33^-1( -_3 β∂_β p  +  ıω v_3 + _3 if_i ),
    ∂_3 v_3   =   ıωκ p - ∂_α v_α + q.


Einstein's summation convention applies also to repeated Greek subscripts (which take on the values 1 and 2).
The particle velocity v_α needs to be eliminated from equation (<ref>).
From equation  (<ref>) we obtain


    v_α=1/ıω(_αβ∂_β p + _α 3∂_3 p - _α if_i ).


Substituting equation (<ref>)  into equation  (<ref>), using equation (<ref>), we obtain


    ∂_3 v_3   =   ıωκ p - 1/ıω∂_α(_αβ∂_β p + _α 3∂_3 p - _α if_i ) + q
       =   ıωκ p - 1/ıω∂_α(_αβ∂_β p + 
    _α 3_33^-1( -_3 β∂_β p  +  ıω v_3 + _3 if_i ) - _α if_i ) + q.


We define


    b_α i=_α i-  _α 3_33^-1_3i,


with b_α 3=0 and b_αβ=b_βα on account of _hi=_ih.
Equations (<ref>) and (<ref>) have the form of equations (<ref>) and (<ref>), with  and 
defined in equation (<ref>), 1 × 1 vectors _1 and _2 defined as


    _1= _33^-1_3if_i,  _2=1/ıω∂_α( b_αβ f_β) + q,


and 1× 1  operator matrices _11, _12, _21 and _22 defined as


    _11   =    -_33^-1_3 β∂_β,
    _12   =   ıω_33^-1, 
    _21   =   ıωκ-1/ıω∂_αb_αβ∂_β, 
    _22   =    -∂_α_α 3_33^-1.


The notation in the right-hand side of these equations  should be understood in the sense that differential operators 
act on all factors to the right of it.  For example, operator ∂_α b_αβ∂_β, applied via equation (<ref>) to p, stands for ∂_α (b_αβ∂_β p), etc.


Operators _11, _12, _21 and _22 
obey the symmetry relations (<ref>) - (<ref>).


We define adjoint acoustic medium parameters as  κ̅= κ^*, _hi = _hi^* and hence b̅_αβ=b_αβ^*.
Operators _11, _12, _21 and _22 in the adjoint medium
are defined as in equations (<ref>) - (<ref>), but with  κ, _hi and b_αβ replaced by  κ̅, _hi  and b̅_αβ, respectively.
These operators obey relations (<ref>) - (<ref>).


For the special case of an isotropic fluid we have _hi=1/ρδ_hi, with ρ denoting the mass density of the isotropic fluid. 
For this situation  equations (<ref>) - (<ref>) reduce to the well-known expressions


    _11   =   _22=0, 
    _12   =   ıωρ, 
    _21   =   ıωκ-1/ıω∂_α1/ρ∂_α,


<cit.>.




§ QUANTUM MECHANICAL WAVES


Schrödinger's wave equation for a particle with mass m in a potential V=V() is given by <cit.>


    ıħ∂_t ψ  = - ħ^2/2m∂_i∂_i ψ + Vψ,


where ψ=ψ(, t) is the wave function and ħ=h/2π, with h Planck's constant.
We use equation (<ref>) to transform this equation to the space-frequency domain, which means we can replace ∂_t by -ıω. 
This gives


    ħωψ  = - ħ^2/2m∂_i∂_i ψ + Vψ,


with ψ=ψ(,ω). To account for losses, we replace V() by a complex-valued, frequency-dependent function V(,ω).

The quantities ψ and  2ħ/mı∂_3ψ constitute the probability current density  in the x_3-direction, via


    =2ħ/mı(ψ^*∂_3ψ-ψ∂_3ψ^*).


We choose these quantities for the 1× 1 vectors  and  in equations (<ref>) and (<ref>), hence


    =ψ,  =2ħ/mı∂_3ψ .


To arrive at a set of equations for these quantities, 
we first recast equation (<ref>) (using the fact that ħ and m are constants) as


    ∂_3(2ħ/mı∂_3ψ) = 4ı(ω-V/ħ)ψ-2ħ/mı∂_α∂_αψ.


This equation, together with the trivial equation 


    ∂_3ψ= mı/2ħ(2ħ/mı∂_3ψ),


have the form of equations (<ref>) and (<ref>), with  and 
defined in equation (<ref>), _1=_2=0, 
and 1× 1  operator matrices _11, _12, _21 and _22 defined as


    _11   =   _22=0, 
    _12   =   mı/2ħ, 
    _21   =   4ı(ω-V/ħ)-2ħ/mı∂_α∂_α.


Operators _12 and _21 
obey the symmetry relations (<ref>) and (<ref>).


We define the adjoint  potential as V̅=V^*. Operators _12 and _21 for the adjoint potential
obey relations (<ref>) and (<ref>).




§ ELECTROMAGNETIC WAVES


In the space-frequency domain, the Maxwell equations for electromagnetic wave propagation read
<cit.>


    -ıω D_i + J_i - ϵ_ijk∂_j  H_k   =    -  J_i^e,
    
    -ıω B_k + ϵ_klm∂_l  E_m   =    -  J_k^m,


where E_m= E_m( x,ω) is the electric field strength, H_k=H_k( x,ω) the magnetic field strength, 
D_i= D_i( x,ω) the electric flux density, B_k=B_k( x,ω) the magnetic flux density, 
J_i= J_i( x,ω) the induced electric current density,
J_i^e= J_i^e( x,ω) and J_k^m= J_k^m( x,ω) are source functions in terms of
external electric and magnetic current densities and, finally, ϵ_ijk is the alternating tensor (or Levi-Civita
tensor), with ϵ_123=ϵ_312=ϵ_231=-ϵ_213=-ϵ_321=-ϵ_132=1 and all other elements being equal to 0.
The constitutive relations for an inhomogeneous, anisotropic, dissipative medium are given by


    D_i    =   ε_ik E_k= ε_0ε_ r,ik E_k,
    
    B_k    =   μ_km H_m = μ_0μ_ r,kmH_m,
    
    J_i    =   σ_ik E_k,


where ε_ik=ε_ik( x,ω), μ_km=μ_km( x,ω) and
σ_ik=σ_ik( x,ω) are the permittivity, permeability and
conductivity tensors, respectively. The subscripts 0 refer to the
parameters in vacuum and the subscripts r denote relative
parameters for the anisotropic medium. 
These tensors obey the symmetry relations ε_ik=ε_ki, μ_km=μ_mk and
σ_ik=σ_ki, respectively.
Substituting the constitutive relations (<ref>) - (<ref>) into Maxwell's
electromagnetic field equations (<ref>) and (<ref>) yields


    -ıω E_ik E_k - ϵ_ijk∂_j  H_k   =    -  J_i^e,
    
    -ıωμ_km H_m + ϵ_klm∂_l  E_m   =    -  J_k^m,


with


    E_ik=ε_ik- σ_ik/ıω.


A matrix-vector wave equation for electromagnetic waves in an isotropic stratified medium is given by <cit.> and <cit.>.
This has been extended for an anisotropic stratified medium by <cit.>. Here we derive the matrix-vector wave equation for electromagnetic waves
in a  3D inhomogeneous, anisotropic, dissipative medium.

The quantities 


    E_0=[ E_1; E_2 ]     H_0=[   H_2; - H_1 ]


constitute the power-flux density  in the x_3-direction, via


    =( E_0^† H_0+ H_0^† E_0)= (E_1^*H_2-E_2^*H_1+H_2^*E_1-H_1^*E_2).


We choose these quantities for the 2× 1 vectors  and  in equations (<ref>) and (<ref>), hence


    = E_0,  = H_0.


To arrive at a set of equations for these quantities, 
we need to eliminate the remaining wave field quantities, E_3 and H_3, from equations (<ref>) and (<ref>).
We start by rewriting these equations as


    -ıω_1   E_0 -ıω_3E_3  + ∂_3  H_0-
     H_3   =   -  J_0^e,
    
    -ıω_3^t  E_0 -ıω E_33 E_3 +^t  H_0
       =   - J_3^e,
    
    -ıω_1   H_0 -ıω_3H_3  + ∂_3  E_0-
     E_3   =   -  J_0^m,
    
    -ıω_3^t  H_0 -ıωμ_33 H_3 +^t  E_0
       =   - J_3^m,


with


    _1=[ E_11 E_12; E_12 E_22 ],  _3=[ E_13; E_23 ],  _1=[  μ_22 -μ_12; -μ_12  μ_11 ],  _3=[  μ_23; -μ_13 ],



    =[ ∂_1; ∂_2 ],  =[  ∂_2; -∂_1 ],   J_0^e=[ J_1^e; J_2^e;       ],   J_0^m=[   J_2^m; - J_1^m;         ].


Equation (<ref>) implies 


    ^t=[  -∂_1 - ∂_2 ],  ^t=[ -∂_2  ∂_1 ].


We separate the derivatives in the x_3-direction from the lateral derivatives in equations (<ref>) and (<ref>), according to


    ∂_3  E_0    =   ıω_1  H_0 +ıω_3H_3+
     E_3-   J_0^m,
    ∂_3  H_0   =   ıω_1   E_0+ıω_3 E_3 +  H_3-  J_0^e.


The field components E_3 and H_3 need to be eliminated. 
From equations (<ref>) and (<ref>) we obtain


    E_3   =    E_33^-1(- _3^t  E_0  +1/ıω^t  H_0+1/ıω J_3^e),
    
    H_3   =   μ_33^-1(- _3^t  H_0  +1/ıω^t  E_0+1/ıω J_3^m).


Substituting equations (<ref>) and (<ref>) into equations (<ref>) and (<ref>) we obtain


    ∂_3  E_0    =   (_3μ_33^-1^t- E_33^-1_3^t )  E_0
    +(ıω_1 - ıω_3μ_33^-1_3^t +1/ıω E_33^-1^t)  H_0
    
           
    +1/ıω( E_33^-1J_3^e)
    - J_0^m+_3μ_33^-1J_3^m,
    ∂_3  H_0   =   (ıω_1 - ıω_3 E_33^-1_3^t +1/ıωμ_33^-1^t)  E_0
    +(_3 E_33^-1^t - μ_33^-1_3^t)  H_0
          -  J_0^e+ _3 E_33^-1J_3^e + 1/ıω (μ_33^-1J_3^m).


Equations (<ref>) and (<ref>) have the form of equations (<ref>) and (<ref>), with  and 
defined in equation (<ref>), 2 × 1 vectors _1 and _2 defined as


    _1   =   1/ıω( E_33^-1J_3^e) -   J_0^m+_3μ_33^-1J_3^m,
    _2   =    -  J_0^e+ _3 E_33^-1J_3^e + 1/ıω (μ_33^-1J_3^m),


and 2× 2  operator matrices _11, _12, _21 and _22 defined as


    _11   =   _3μ_33^-1^t- E_33^-1_3^t,
    _12   =   ıω(_1 - _3μ_33^-1_3^t) +1/ıω E_33^-1^t,
    _21   =   ıω( _1 - _3 E_33^-1_3^t) +1/ıωμ_33^-1^t,
    _22   =   _3 E_33^-1^t - μ_33^-1_3^t.


These operators obey the symmetry relations (<ref>) - (<ref>).

We define adjoint electromagnetic medium parameters as ε̅_ik=ε_ik^*,
μ̅_km=μ_km^* and σ̅_ik=-σ_ik^*.
Using equation (<ref>) it follows that E̅_ik= E_ik^*.
Similar relations hold for _1, _3, _1 and _3, which contain the parameters E_ik and μ_km. 
Operators _11, _12, _21 and _22 in the adjoint medium
obey relations (<ref>) - (<ref>).

Explicit expressions for the operator matrices in an isotropic medium are given in the supplemental material, Appendix I.




§ ELASTODYNAMIC WAVES



In the space-frequency domain, the elastodynamic equations of motion and deformation read  <cit.>


    -ıω_i - ∂_jτ_ij   =    f_i,
    ıω e_kl+ (∂_k v_l+∂_l v_k)    =    h_kl,








where _i=_i(,ω) is the momentum density,

τ_ij=τ_ij(,ω) the stress tensor,
e_kl=e_kl(,ω) the strain tensor,
v_k=v_k(,ω) the particle velocity,



and f_i=f_i(,ω) and h_kl=h_kl(,ω) are source functions in terms of external force density and deformation-rate density, respectively. 
The stress, strain and  deformation-rate tensors obey the  symmetry relations  τ_ij=τ_ji, e_kl=e_lk and h_kl=h_lk.
The constitutive relations for an inhomogeneous, anisotropic, dissipative solid are given by


    _i    =   ρ_ijv_j,
    
    e_kl   =   s_klmnτ_mn,


where ρ_ij=ρ_ij(,ω) and  s_klmn=s_klmn(,ω) are the mass density and compliance tensors, respectively.
These tensors obey the symmetry relations  ρ_ij=ρ_ji and
s_klmn=s_lkmn=s_klnm=s_mnkl, respectively <cit.>.
Substituting the constitutive relations (<ref>) and  (<ref>) into equations (<ref>) and (<ref>) yields 


    -ıωρ_ijv_j - ∂_jτ_ij   =    f_i,
    ıω s_klmnτ_mn+ (∂_k v_l+∂_l v_k)    =   h_kl.




We introduce the stiffness tensor
c_ijkl=c_ijkl( x,ω) as the inverse of the compliance tensor s_klmn, according to


    c_ijkls_klmn=s_ijklc_klmn=
    (δ_imδ_jn+δ_inδ_jm).


The stiffness tensor obeys the symmetry relation 
c_ijkl=c_jikl=c_ijlk=c_klij.
Multiplying all terms in equation (<ref>) by c_ijkl, using the symmetry relations
τ_ij=τ_ji and  c_ijkl=c_ijlk, we obtain an alternative form of 

equation (<ref>), according to


    ıωτ_ij+ c_ijkl∂_l v_k   =   c_ijklh_kl.


A matrix-vector wave equation for elastodynamic  waves in an inhomogeneous anisotropic medium is given by <cit.>.
Here we review this derivation, which also serves as a starting point for the
derivation of the matrix-vector wave equations for  poroelastodynamic waves (Appendix <ref>), piezoelectric waves (Appendix <ref>) and 
seismoelectric waves (Appendix <ref>).

The quantities -_3 and v (which are 3× 1 vectors, with (_3)_i=τ_i3 and ( v)_i=v_i)
constitute the power-flux density  in the x_3-direction, via


    =(-_3^† v- v^†_3)=  (-τ_i3^*v_i-v_i^*τ_i3).


We choose these quantities for the 3× 1 vectors  and  in equations (<ref>) and (<ref>), hence


    =-_3,  = v.


To arrive at a set of equations for these quantities, 
we need to eliminate the remaining wave field quantities, 3× 1 vectors _1 and _2 (with (_j)_i=τ_ij), from equations (<ref>) and (<ref>).
We start by rewriting these equations as


    -ıω v - ∂_j_j    =    f,
    ıω_j+_jl∂_l  v   =   _jl h_l,


where  and _jl are 3× 3 matrices, with  ()_ij=ρ_ij, ^t=, (_jl)_ik=c_ijkl,  _jl^t=_lj,
and where f and h_l are 3× 1 vectors, with ( f)_i=f_i and ( h_l)_k=h_kl.
We separate the derivatives in the x_3-direction from the lateral derivatives in equations (<ref>) and (<ref>), according to


    -∂_3_3    =   ıω v + ∂_α_α+  f,
    ∂_3  v   =   _33^-1(-ıω_3 -_3β∂_β v +_3l h_l).


The field components _1 and _2 need to be eliminated. 
From equation (<ref>) we obtain


    _α=-1/ıω(_αβ∂_β v +_α 3∂_3  v - _α l h_l).


Substituting equation (<ref>) into (<ref>) and the result into equation  (<ref>), we obtain


    -∂_3_3    =   ∂_α(_α 3_33^-1_3)
    +ıω v-1/ıω∂_α(_αβ∂_β v -_α l h_l) +  f,


with


    _α l=_α l-
    _α 3_33^-1_3l,


where _α 3= O and _αβ^t=_βα on account of _jl^t= _lj.
Equations  (<ref>) and (<ref>) have the form of equations (<ref>) and (<ref>), with  and 
defined in equation (<ref>), 3 × 1 vectors _1 and _2 defined as


    _1   =    f+1/ıω∂_α(_αβ h_β),
    _2   =   _33^-1_3l h_l,


and 3× 3  operator matrices _11, _12, _21 and _22 defined as


    _11   =    -∂_α_α 3_33^-1,
    _12   =   ıω-1/ıω∂_α_αβ∂_β,
    _21   =   ıω_33^-1,
    _22   =   -_33^-1_3β∂_β.


These operators obey the symmetry relations (<ref>) - (<ref>).

We define adjoint elastodynamic medium parameters as c̅_ijkl=c_ijkl^* and ρ̅_ij=ρ_ij^*.
Similar relations hold for _jl, _αβ and , which contain the parameters c_ijkl and ρ_ij. 
Operators _11, _12, _21 and _22 in the adjoint medium
obey relations (<ref>) - (<ref>).

Explicit expressions for the operator matrices in an isotropic medium are given in the supplemental material, Appendix J.



§ POROELASTODYNAMIC WAVES


In the space-frequency domain, the basic equations for poroelastodynamic wave propagation in an inhomogeneous, anisotropic, dissipative,
fluid-saturated porous solid read  <cit.>


    -ıωρ_ij^bv_j^s -ıωρ_ij^fw_j - ∂_jτ_ij^b   =    f_i^b,
    
    -ıω/η k_ijρ_jl^fv_l^s +w_i + 1/η k_ij∂_jp^f   =   1/η k_ij f_j^f,
    ıωτ_ij^b+ _̧ijkl∂_lv_k^s +  C_ij∂_kw_k   =   _̧ijkl_kl+C_ij,
    
    -ıωp^f+ C_kl∂_lv_k^s+ M∂_kw_k   =   C_kl_kl+M,


with


    w_j   =   ϕ( v_j^f- v_j^s),
    
    v_j^b   =   ϕ v_j^f+(1-ϕ)v_j^s=v_j^s+w_j,
    τ_ij^b   =   ϕτ_ij^f+(1-ϕ)τ_ij^s=-ϕδ_ijp^f+(1-ϕ)τ_ij^s,
    
    f_i^b   =   ϕ f_i^f+(1-ϕ)f_i^s,
    ρ_ij^b   =   ϕρ_ij^f+(1-ϕ)ρ_ij^s.


Superscripts b, f and s stand for bulk, fluid and solid, respectively. The wave field quantity
v_j= v_j( x,ω) is the averaged particle velocity in the bulk, fluid or solid (depending on the superscript),
w_j= w_j( x,ω) is the filtration velocity, τ_ij=τ_ij( x,ω)  the averaged stress in the bulk, fluid or solid
and p^f= p^f( x,ω) the averaged fluid pressure.
The stress tensors are symmetric, i.e., τ_ij=τ_ji.
The medium parameter ρ_ij=ρ_ij( x,ω)
is the mass density of the bulk, fluid or solid (depending on the superscript). 
Furthermore,  k_ij=k_ij( x,ω) is the dynamic permeability tensor,
η=η( x,ω) is the fluid viscosity parameter and ϕ=ϕ( x) the porosity.
Moreover, _̧ijkl= _̧ijkl( x,ω), C_ij= C_ij( x,ω) and M= M( x,ω) are
stiffness parameters of the porous solid. The medium parameters obey the following symmetry relations
ρ_ij=ρ_ji, k_ij=k_ji, _̧ijkl= _̧jikl= _̧ijlk= _̧klij and  C_ij= C_ji.
The source function f_i= f_i( x,ω) is the volume density of external force on the bulk, fluid or solid.
For many source types the forces on the bulk and fluid are equal but in the following they will be treated distinctly.
The source functions _kl=_kl( x,ω) and =( x,ω) are the volume densities of external deformation rate on the bulk and volume-injection rate in the fluid
<cit.>. The deformation rate tensor is symmetric, i.e., _kl=_lk.
For later convenience, we eliminate ∂_kw_k from equation
(<ref>), using equation (<ref>). This yields


    ıωτ_ij^b+ _ijkl∂_lv_k^s +ıω/ MC_ijp^f=_ijkl_kl,


with  _ijkl= _ijkl( x,ω) defined as 


    _ijkl= _̧ijkl-1/ MC_ijC_kl.


A matrix-vector wave equation for normal-incidence poroelastodynamic waves in a stratified isotropic medium is given by <cit.> and <cit.>.
This has been extended for oblique-incidence poroelastodynamic waves in a stratified anisotropic medium, separately for P-SV and SH propagation, by <cit.>.
Here we derive the matrix-vector wave equation for poroelastodynamic waves in a 3D inhomogeneous, anisotropic, dissipative, fluid-saturated porous solid. 

The quantities -_3^b, p^f, v^s and w_3 (with (_j^b)_i=τ_ij^b and ( v^s)_i=v_i^s)
constitute the power-flux density  in the x_3-direction, via


    =   (-(_3^b)^† v^s+ p^f*w_3 - ( v^s)^†_3^b + w_3^*p^f)
       =   (-τ_i3^b*v_i^s + p^f*w_3 - v_i^s*τ_i3^b + w_3^*p^f).


We choose these quantities for the 4× 1 vectors  and  in equations (<ref>) and (<ref>), hence


    =[ - _3^b;    p^f ],
      =[ v^s; w_3 ].


To arrive at a set of equations for these quantities, 
we need to eliminate the remaining wave field quantities  from equations  (<ref>), (<ref>), (<ref>) and (<ref>).
We start by rewriting these equations as


    -ıω^b v^s -ıω^f_j w_j -∂_j_j^b    =    f^b,
    
    -ıω/η^f v^s +_j w_j + 1/η_j∂_jp^f   =   1/η_j f_j^f,
    ıω_j^b+  _jl∂_l  v^s+ ıω/ M_jp^f   =   _jl_l,
    
    -ıωp^f+ _l^t ∂_l v^s+  M ∂_kw_k    =   _l^t_l +M,


where ,  and _jl are 3× 3 matrices, with  ()_ij=ρ_ij, ^t=,  ()_ij=k_ij, ^t=, (_jl)_ik=_ijkl,  
_jl^t=_lj, and where _j, f^b, _l and _j are 3× 1 vectors, with (_j)_i=C_ij, ( f^b)_i=f_i^b, (_l)_k=_kl, 
and (_j)_i=δ_ij.

Equations (<ref>) - (<ref>) form the starting point for deriving matrix-vector equations in the form of
equations (<ref>) and (<ref>), with  and  defined in equation  (<ref>).
The other quantities (_1^b, _2^b, w_1 and w_2) need to be eliminated. 
The detailed derivation can be found in the supplemental material, Appendix K. 
The 4× 1 vectors _1 and _2 are defined as


    _1     =   [ ıω/η^f_α_α^t( I+1/_3_3^t^-1_γ_γ^t)_jf_j^f
    -ıω/η^f_α_α^t_3f_3^f
    + f^b+1/ıω∂_α(_αβ_β);                                                     1/(-_3^t^-1_α_α^t_jf_j^f+f_3^f) ],
    _2    =   [                                                             _33^-1_3l_l; -∂_α(1/η_α^t( I+
    1/_3_3^t^-1_γ_γ^t)_jf_j^f-1/η_α^t_3f_3^f )+
    1/M_α^t_α+ ]


and the 4× 4 operator matrices _11, _12, _21 and _22 as


    _11   =   [ _11^11 _11^12; _11^21 _11^22 ],  _12 = [     _12^11     _12^12; (_12^12)^t     _12^22 ],
    _21   =   [     _21^11     _21^12; (_21^12)^t     _21^22 ],  _22 =-_11^t.


Here


    _11^11   =    -∂_α_α 3_33^-1,
    _11^12   =   -ıω/η^f( I-1/_3_3^t^-1_3_3^t)_β∂_β 
    -∂_α1/M_α,
    _11^21   =    0^t,
    _11^22   =   -1/_3^t^-1_3_3^t_β∂_β,


    _12^11   =   ıω^b-1/ıω∂_α_αβ∂_β
    -ω^2/η^f_α_α^t(-1/_3_3^t^-1_3_3^t)_β_β^t^f,
    _12^12   =   ıω/^f_3_3^t^-1_3,
    _12^22   =   -η/_3^t^-1_3,
    _21^11   =   ıω_33^-1,
    _21^12   =   -ıω/M_33^-1_3,
    _21^22   =   ıω/M^2_3^t_33^-1_3 +ıω/M
    +∂_α1/η(_α^t_β - 1/_α^t_3_3^t^-1_3_3^t_β)∂_β,


with 0 denoting a  zero vector and


    _αβ   =   _αβ-
    _α 3_33^-1_3β,
    _α   =   _α-  _α 3_33^-1_3,
    =   1-_3^t^-1_α_α^t_3.


Operators _11, _12, _21 and _22 
obey the symmetry relations formulated in equations (<ref>) - (<ref>).

We defined the adjoints of the medium parameters  c_ijkl and ρ_ij in Appendix <ref>
(where ρ_ij now has superscript b or f). Moreover, we define
k̅_ij=-k_ij^*, η̅= η^*, C̅_ij=C_ij^* and M̅= M^*.
Similar relations hold for  _jl, ^b, ^f,  and _j, 
which contain the parameters _ijkl= _̧ijkl-1/ MC_ijC_kl, ρ_ij^b,  ρ_ij^f, k_ij and C_ij.
Operators _11, _12, _21 and _22 in the adjoint medium
obey relations (<ref>) - (<ref>).

Explicit expressions for the operator matrices in an isotropic medium are given in the supplemental material, Appendix K.



§ PIEZOELECTRIC WAVES


Piezoelectric waves are governed by the equations for electromagnetic waves (Appendix <ref>) and elastodynamic waves (Appendix <ref>), in which two of the constitutive 
relations need to be modified to account for the coupling between the two wave types. For piezoelectric waves, the modified constitutive relations are <cit.>


    D_i    =   ε_ik E_k + d_ijkτ_jk,
    
    e_kl   =   d_klmE_m+s_klmnτ_mn.


The field quantities and medium parameters (except d_ijk) have been defined in Appendices <ref> and <ref>.
Parameters ε_ik  in equation (<ref>) and s_klmn in equation (<ref>)
are defined under constant stress and constant electric field, respectively.
The coupling tensor d_ijk=d_ijk(,ω) obeys the symmetry relation d_ijk=d_jik=d_ikj.
Equation (<ref>) replaces constitutive relation (<ref>) and is substituted, together with constitutive relations (<ref>) and (<ref>),
into the Maxwell equations (<ref>) and (<ref>).
Equation (<ref>) replaces stress-strain relation (<ref>)
and is substituted, together with constitutive relation (<ref>), into equations (<ref>) and (<ref>).  
Subsequently, all terms in the latter equation are multiplied by c_ijkl, using equation (<ref>) as well as the symmetry relations
τ_ij=τ_ji and  c_ijkl=c_ijlk.


The basic equations for coupled electromagnetic and elastodynamic waves thus read


    -ıω E_ikE_k -ϵ_ijk∂_j H_k-ıω d_ijkτ_jk   =   -J_i^e,
    
    -ıωμ_kmH_m+ϵ_klm∂_l E_m   =   -J_k^m,
    
    -ıωρ_ijv_j -∂_jτ_ij    =    f_i,
    ıωτ_ij+ c_ijkl∂_l v_k+ ıω c_ijkld_klmE_m   =    c_ijklh_kl,


with E_ik=ε_ik- σ_ik/ıω.

A matrix-vector wave equation in the quasi-static approximation for  2D piezoelectric waves  in an anisotropic stratified medium  is given 
by <cit.>, <cit.> and <cit.>. 
Here we derive the exact matrix-vector wave equation for piezoelectric waves in a 3D inhomogeneous, anisotropic, dissipative,
piezoelectric medium. 

The quantities -_3, E_0, v and H_0
(with E_0 and H_0 defined in Appendix <ref> and _3  and v defined in Appendix <ref>)
constitute the power-flux density  in the x_3-direction, via


    =   (-_3^† v+ E_0^† H_0- v^†_3+ H_0^† E_0)
       =   (-τ_i3^*v_i+E_1^*H_2-E_2^*H_1-v_i^*τ_i3+H_2^*E_1-H_1^*E_2).


We choose these quantities for the 5× 1 vectors  and  in equations (<ref>) and (<ref>), hence


    =[ - _3;  E_0 ],
      =[   v; H_0 ].


To arrive at a set of equations for these quantities, 
we need to eliminate the remaining wave field quantities  from equations (<ref>) - (<ref>).
Using the notation introduced in Appendices <ref> and <ref>,
we rewrite equations (<ref>) - (<ref>) as


    -ıω_1   E_0 -ıω_3E_3  + ∂_3  H_0- H_3-ıω D_1k^t_k   =   -  J_0^e,
    
    -ıω_3^t  E_0 -ıω E_33 E_3 +^t  H_0-ıω D_3k^t_k   =   - J_3^e,
    
    -ıω_1   H_0 -ıω_3H_3  + ∂_3  E_0- E_3   =   -  J_0^m,
    
    -ıω_3^t  H_0 -ıωμ_33 H_3 +^t  E_0   =   - J_3^m,
    
    -ıω v - ∂_j_j    =    f,
    ıω_j+_jl∂_l  v+ıω_jl( D_1l E_0+ D_3lE_3)   =   _jl h_l,


with


    D_1k=[ d_11k d_12k; d_21k d_22k; d_31k d_32k ],   D_3k=[ d_13k; d_23k; d_33k ].


Equations (<ref>) - (<ref>) form the starting point for deriving 
matrix-vector equations in the form of equations (<ref>) and (<ref>), with  and  defined in equation (<ref>).
The other quantities (_1,  _2, E_3 and H_3) need to be eliminated.
The detailed derivation can be found in the supplemental material, Appendix L. 
The 5 × 1 vectors _1 and _2 are defined as


    _1   =   [       f+1/ıω∂_α( _αβ' h_β)  - 1/ıω∂_α(( E_33')^-1_αβ D_3β J_3^e); 1/ıω(( E_33')^-1(J_3^e- D_3α^t_αβ h_β))-   J_0^m+_3μ_33^-1 J_3^m;                                                                  ],
    _2   =   [                            _33^-1(_3m' h_m-( E_33')^-1_3l D_3l J_3^e); -  J_0^e  + ( E_33')^-1_3'J_3^e+1/ıω(μ_33^-1J_3^m)+ ( D_1α')^t_αβ h_β ]


and the 5× 5  operator matrices _11, _12, _21 and _22 
having the form defined in equations (<ref>) and (<ref>), where 


    _11^11   =   -∂_α(_3α')^t_33^-1,
    _11^12   =   -∂_α_αβ D_1β',
    _11^21   =   ( E_33')^-1 D_3k^t_k 3_33^-1,
    _11^22   =   _3μ_33^-1^t - ( E_33')^-1(_3')^t,
    _12^11   =   ıω-1/ıω∂_α_αβ' ∂_β, 
    _12^12   =   -1/ıω∂_α( E_33')^-1_αβ D_3β^t,
    _12^22   =   ıω(_1- _3μ_33^-1_3^t)+1/ıω( E_33')^-1^t,
    _21^11   =   ıω(_33^-1 - ( E_33')^-1_33^-1_3l D_3l D_3k^t_k 3_33^-1),
    _21^12   =   -ıω_33^-1_3l D_1l',


    _21^22   =   ıω(_1' -_3'( E_33')^-1(_3')^t)+1/ıωμ_33^-1^t,


with


    _αβ   =   _αβ-_α 3_33^-1_3β,
    _1'   =   _1- D_1α^t_αβ D_1β,
    _3'   =   _3- D_1α^t_αβ D_3β,
     E_33'   =    E_33- D_3α^t_αβ D_3β,
    _αβ'   =   _αβ+( E_33')^-1_αγ D_3γ D_3δ^t _δβ,
    _3m'   =   _3m+ ( E_33')^-1_3l D_3l D_3α^t_α m,
     D_1l'    =    D_1l  - ( E_33')^-1 D_3l (_3')^t,


with _α 3= O. Operators _11, _12, _21 and _22 
obey the symmetry relations formulated in equations (<ref>) - (<ref>).

We defined the adjoints of the medium parameters E_ik,  μ_km, c_ijkl and ρ_ij in Appendices <ref> and <ref>.
Moreover, we define d̅_ijk=d_ijk^*.

Similar relations hold for _1, _3, _1, _3,  _jl, , D_1k and D_3k, 
which contain the parameters E_ik, μ_km, c_ijkl, ρ_ij and d_ijk.
Operators _11, _12, _21 and _22 in the adjoint medium obey relations (<ref>) - (<ref>).



§ SEISMOELECTRIC WAVES


Seismoelectric waves are governed by the equations for electromagnetic waves (Appendix <ref>) and poroelastodynamic waves (Appendix <ref>), in which two of 
the constitutive relations need to be modified to account for the coupling  between the two wave types. In this appendix we consider an isotropic medium
(the derivation for the anisotropic situation is disproportionally long). For seismoelectric waves, the modified constitutive relations are  <cit.>


    J_i   =   σ E_i + L(-∂_ip^f+ıωρ^fv_i^s+f_i^f),
    
    w_i   =   L E_i + k/η(-∂_ip^f+ıωρ^fv_i^s+f_i^f).


The field quantities, sources and medium parameters (except L)  have been defined in Appendices
<ref> and <ref> (except that tensors are now replaced by scalars).
Here  L=L( x,ω) accounts for the coupling  between the elastodynamic and
electromagnetic waves and vice versa.
Equations (<ref>) and (<ref>) contain the same  coupling coefficient L  (due to Onsager's reciprocity relation, <cit.>). 
Equation (<ref>) replaces the isotropic version of 
constitutive relation (<ref>) and is substituted, together with the isotropic versions of constitutive relations (<ref>) and (<ref>),
into the Maxwell equations (<ref>) and (<ref>). Equation (<ref>) replaces the isotropic version of equation (<ref>).
The basic equations for coupled electromagnetic and poroelastodynamic waves thus read


    -ıωρ^bv_i^s -ıωρ^fw_i - ∂_jτ_ij^b   =    f_i^b,
    
    -ıωρ^fv_i^s +η/k(w_i - LE_i) + ∂_ip^f   =     f_i^f,
    ıωτ_ij^b+ _̧ijkl∂_lv_k^s +  Cδ_ij∂_kw_k   =   _̧ijkl_kl+Cδ_ij,
    
    -ıωp^f+Cδ_kl∂_lv_k^s+ M∂_kw_k   =   Cδ_kl_kl+M,
    
    -ıωε E_i+σ E_i +L(-∂_ip^f+ıωρ^fv_i^s) - ϵ_ijk∂_j H_k   =   - J_i^e - Lf_i^f,
    
    -ıωμ H_k + ϵ_klm∂_l E_m   =    - J_k^m.


For the isotropic medium we have


    _̧ijkl   =   ( K_G-2/3 G_ fr)δ_ijδ_kl+ G_ fr(δ_ikδ_jl+δ_ilδ_jk),


where G_ fr is the
shear modulus of the framework of the grains when the fluid is absent and K_G is the Gassmann modulus <cit.>. 
The permittivity and permeability are defined as ε=ε_0ε_ r and μ=μ_0μ_ r. 
The subscripts 0 refer to the parameters in vacuum and the subscripts r denote
relative parameters. For ε_ r and μ_ r we have  <cit.> 


    ε_ r   =   ϕ/α_∞(κ^f-κ^s)+κ^s,
    μ_ r   ≈   1,


where κ^f and κ^s are the dielectric parameters of the fluid and solid, respectively, and α_∞ is the tortuosity at infinite frequency.

For later convenience, we eliminate ∂_kw_k from equation
(<ref>), using equation (<ref>). This yields


    ıωτ_ij^b+ _ijkl∂_lv_k^s +ıω C/ Mδ_ijp^f= _ijkl_kl,


with _ijkl= _ijkl( x,ω) defined  as


    _ijkl= _̧ijkl- C^2/ Mδ_ijδ_kl.


Also for later convenience, we add L times equation (<ref>) to
equation (<ref>) in order to compensate for the term L(-∂_ip^f+ıωρ^fv_i^s). This yields


    -ıω E  E_i +η/kL w_i- ϵ_ijk∂_j H_k   =   - J_i^e,


with 


    E= ε-1/ıω(σ -η/k L^2).



A matrix-vector wave equation for oblique-incidence seismoelectric waves in a stratified isotropic medium, separately for  P-SV-TM and SH-TE propagation, is given by
<cit.>,  <cit.> and <cit.>.
Here we derive the matrix-vector wave equation for a 3D inhomogeneous, isotropic, dissipative, fluid-saturated porous solid.

The quantities -_3^b, p^f, E_0, v^s, w_3 and H_0 
constitute the power-flux density  in the x_3-direction, via


    =   (-(_3^b)^† v^s+ p^f*w_3 + E_0^† H_0 - ( v^s)^†_3^b + w_3^*p^f + H_0^† E_0)
       =   (-τ_i3^b*v_i^s + p^f*w_3 + E_1^*H_2-E_2^*H_1 - v_i^s*τ_i3^b + w_3^*p^f +H_2^*E_1-H_1^*E_2).


We choose these quantities for the 6× 1 vectors  and  in equations (<ref>) and (<ref>), hence


    =[ - _3^b;    p^f;    E_0 ],
      =[ v^s; w_3; H_0 ].


To arrive at a set of equations for these quantities, we need to eliminate the remaining wave field quantities from
equations (<ref>), (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>).
We start by rewriting these equations  as


    -ıωρ^b v^s -ıωρ^f_j w_j -∂_j_j^b    =    f^b,
    
    -ıωρ^f_i^t v^s +η/k( w_i -L(_i^t E_0+δ_3i E_3))+ ∂_ip^f   =   f_i^f,
    ıω_j^b+  _jl∂_l  v^s+ıω C/ M_jp^f   =   _jl_l,
    
    -ıωp^f+  C_l^t ∂_l v^s+ M ∂_kw_k    =   C_l^t _l+M,
    
    -ıω E E_0 +η/ kL_α w_α+∂_3  H_0- H_3   =   -  J_0^e,
    
    -ıω E  E_3 +η/ kL w_3+^t  H_0   =   - J_3^e,
    
    -ıωμ H_0 + ∂_3  E_0 - E_3   =   - J_0^m,
    
    -ıωμ H_3 +^t  E_0   =    - J_3^m,


with most of the vectors and matrices defined in Appendices <ref> and  <ref>. In addition, _i is a 2× 1 unit vector, with (_i)_β=δ_β i.
Equations (<ref>) - (<ref>) form the starting point for deriving matrix-vector equations in the form of
equations (<ref>) and (<ref>), with  and  defined in equation  (<ref>).
The other quantities (_1^b, _2^b, E_3, w_1, w_2 and H_3)  need to be eliminated. 
The detailed derivation can be found in the supplemental material, Appendix M. 
The 6× 1 vectors _1 and _2 are defined as


    _1     =    [ f^b+ıωρ^fk/η_α f_α^f+1/ıω∂_α(_αβ_β);             1/ıω Eη/ kLJ_3^e+ f_3^f;               - J_0^m+(1/ıω EJ_3^e) ],


    _2     =    [                      _33^-1_3l_l;      +1/M_α^t_α-∂_β(  k/η f_β^f); - J_0^e +(1/ıωμ J_3^m)-L_α f_α^f ]


and the 6× 6  operator matrices _11, _12, _21 and _22 as


    _11   =   [ _11^11 _11^12 _11^13;    0^t      0    0^t;      O      0      O ],  _12 = [     _12^11     _12^12          O; (_12^12)^t     _12^22     _12^23;          O (_12^23)^t     _12^33 ],
    _21   =   [     _21^11     _21^12          O; (_21^12)^t     _21^22     _21^23;          O (_21^23)^t     _21^33 ],  _22 = -_11^t.


Here


    _11^11   =    -∂_α_α 3_33^-1,
    _11^12   =   -ıωρ^fk/η_α∂_α-∂_α1/ M_α,
    _11^13   =   ıωρ^fL_α_α^t,
    _12^11   =   -1/ıω∂_α_αβ∂_β+ıω( ρ^b I_3+ıω( ρ^f)^2k/η_α_α^t),
    _12^12   =   ıωρ^f_3,
    _12^22   =   -η/ k(1-1/ıω Eη/ kL^2),
    _12^23   =   1/ıω Eη/ kL^t,
    _12^33   =   ıωμ I_2 +1/ıω E^t,
    _21^11   =   ıω_33^-1,
    _21^12   =   -ıω C/ M_33^-1_3,
    _21^22   =   ıω C^2/ M^2_3^t _33^-1_3 +ıω/ M+∂_β k/η∂_β,
    _21^23   =   -∂_βL_β^t,
    _21^33   =   (ıω E-η/ kL^2) I_2+1/ıωμ^t,


where I_3 is a 3× 3 identity matrix, I_2  a 2× 2 identity matrix and


    _αβ   =   _αβ- _α 3_33^-1_3β,
    _α   =   C(_α-  _α3_33^-1_3).


Operators _11, _12, _21 and _22 
obey the symmetry relations formulated in equations (<ref>) - (<ref>).

We define the adjoints of the medium parameters  ϵ, μ, σ, c_ijkl, ρ (with superscript b or f), k, η, C and M 
similar as in Appendices <ref>, <ref> and <ref> but for the isotropic situation.
Moreover, we define L̅ = -L^*.
Operators _11, _12, _21 and _22 in the adjoint medium
obey relations (<ref>) - (<ref>).

Explicit expressions for the operator matrices are given in the supplemental material, Appendix M.



§ SYMMETRY PROPERTIES OF THE OPERATOR MATRIX IN THE WAVENUMBER-FREQUENCY DOMAIN


We derive symmetry properties of the operator matrix for the special case of a laterally invariant medium (or potential) in the wavenumber-frequency domain.
We define the spatial Fourier transform of a space- and frequency-dependent quantity h(,ω) as


    h̃(k_α,x_3,ω)=∫_ h(,ω) exp(-ı k_α x_α) d^2,


with k_α for α=1,2 representing the horizontal wavenumbers.
Lateral derivatives ∂_α h(,ω) in the space-frequency domain are replaced by products
ı k_αh̃(k_α,x_3,ω) in the wavenumber-frequency domain. We denote the Fourier transform
of ∂_α as ∂_α⇒ı k_α.
Similarly, for an operator matrix  in a laterally invariant medium, containing the differential operator ∂_α, we denote the Fourier transform as 
(∂_α) ⇒(k_α).
Using equations (<ref>) and (<ref>), we find


    {(∂_α)}^t    ⇒   {(-k_α)}^t,
    {(∂_α)}^*    ⇒   {(-k_α)}^*,
    {(∂_α)}^†   ⇒   {(k_α)}^†.


We use equation (<ref>) to transform equation  (<ref>)  to the wavenumber-frequency domain, according to 
∂_3=+.
We find the symmetry properties of (k_α,x_3,ω) by applying equations (<ref>) - (<ref>) to
 the left-hand sides of equations (<ref>) - (<ref>). This gives


    {(-k_α,x_3,ω)}^t   =   - (k_α,x_3,ω),
    {(-k_α,x_3,ω)}^*   =      (k_α,x_3,ω),
    {(k_α,x_3,ω)}^†   =   -    (k_α,x_3,ω).


Unified matrix-vector wave equation, reciprocity and representations: Supplemental material





§ ELECTROMAGNETIC WAVES


The electromagnetic matrix-vector wave equation for an inhomogeneous, anisotropic, dissipative medium is derived in Appendix C in the main paper.
For the special case of an isotropic medium we have E_ik= Eδ_ik and μ_km=μδ_km,  or


    _1=[ E 0; 0 E ],  _3=[ 0; 0 ],   E_33= E,  _1=[ μ 0; 0 μ ],  _3=[ 0; 0 ],  μ_33=μ.


For this situation the wave and source vectors read 


    = [ E_0; H_0 ],    = [ 1/ıω( E^-1J_3^e) -   J_0^m; -  J_0^e+ 1/ıω (μ^-1J_3^m) ],


and the operator matrices reduce to


    _12   =   ıω_1  +1/ıω E^-1^t=
    [ ıωμ-1/ıω∂_1
    1/ E∂_1    -1/ıω∂_1
    1/ E∂_2;    -1/ıω∂_2
    1/ E∂_1 ıωμ-1/ıω∂_2
    1/ E∂_2 ],
    _21   =   ıω_1 +1/ıωμ^-1^t=
    [ ıω E
    -1/ıω∂_2
    1/μ∂_2       1/ıω∂_2
    1/μ∂_1;       1/ıω∂_1
    1/μ∂_2 ıω E -1/ıω∂_1 1/μ∂_1 ]


and _11=_22=, where  is a 2× 2 zero matrix.




§ ELASTODYNAMIC WAVES


The elastodynamic matrix-vector wave equation for an inhomogeneous, anisotropic, dissipative solid is derived in Appendix D in the main paper.
For the special case of an isotropic medium we have


    c_ijkl   =   λδ_ijδ_kl+ μ(δ_ikδ_jl+δ_ilδ_jk),
    ρ_ij   =   ρδ_ij,


with Lamé parameters λ=λ( x,ω) and μ=μ( x,ω) and mass density  ρ=ρ(,ω).
Hence,   the mass density matrix reduces to =ρ I. For the stiffness matrices _jl we have (_jl)_ik=c_ijkl=λδ_ijδ_kl+ μ(δ_ikδ_jl+δ_ilδ_jk), hence


    [  _11=[ K_c   0   0;   0   μ   0;   0   0   μ ],                   _12=
    [ 0 λ 0; μ 0 0; 0 0 0 ],                   _13=
    [ 0 0 λ; 0 0 0; μ 0 0 ],;                                                ;                   _21=
    [ 0 μ 0; λ 0 0; 0 0 0 ], _22=
    [   μ   0   0;   0 K_c   0;   0   0   μ ],                   _23=
    [ 0 0 0; 0 0 λ; 0 μ 0 ],                                                 ]


    [                   _31=
    [ 0 0 μ; 0 0 0; λ 0 0 ],                   _32=
    [ 0 0 0; 0 0 μ; 0 λ 0 ], _33=
    [   μ   0   0;   0   μ   0;   0   0 K_c ],;                                                 ]


with K_c=K_c( x,ω)=λ( x,ω) +2μ( x,ω).

For this situation the operator matrices reduce to

    _11   =   [           0           0 -∂_1
    λ/λ+2μ;           0           0 -∂_2
    λ/λ+2μ;        -∂_1        -∂_2           0 ],
          
    _12   =   [ ıωρ-1/ıω(
    ∂_1ν_1∂_1+
    ∂_2μ∂_2)   -1/ıω(∂_2 μ∂_1+ ∂_1 ν_2∂_2)                             0; -1/ıω(∂_2
    ν_2
    ∂_1+ ∂_1 μ∂_2 ) ıωρ-1/ıω(
    ∂_1μ∂_1+ ∂_2ν_1∂_2)                             0;                             0                             0                           ıωρ ],
          
    _21   =   [    ıω/μ       0       0;       0    ıω/μ       0;       0       0 ıω/λ+2μ ],
          
    _22   =   [          0          0       -∂_1;          0          0       -∂_2; -λ/λ+2μ∂_1 -λ/λ+2μ∂_2          0 ],


where


    ν_1   =   ν_1( x,ω)=
    4μ(λ+μ/λ+2μ),
    ν_2   =   ν_2( x,ω)=
    2μ(λ/λ+2μ).





§ POROELASTODYNAMIC WAVES



Equations (E.14) - (E.17) in Appendix E in the main paper form the starting point for deriving a poroelastodynamic
matrix-vector wave equation 
for the quantities _3^b, p^f,   v^s and w_3 in vector . 
Pre-multiplying all terms in  equation (E.15) by η_i^t^-1, with 3× 1 unit vector (_i)_j=δ_ij, gives


    -ıω_i^t^f v^s +η_i^t^-1_j w_j + ∂_ip^f= f_i^f.


We separate the derivatives in the x_3-direction from the lateral
derivatives in equations (E.14),

(<ref>),

(E.16) and (E.17), according to


    -∂_3_3^b    =   ıω^b v^s
    +ıω^f (_α w_α+_3 w_3 )
    +∂_α_α^b+
     f^b,
    
    ∂_3  p^f    =    ıω_3^t^f v^s -η_3^t^-1(_α w_α+_3 w_3 )+ f_3^f,
    ∂_3  v^s   =   _33^-1(-ıω_3^b-ıω/ M_3 p^f - _3β∂_β v^s
    +_3l_l),
    ∂_3  w_3    =   ıω/ M p^f -  1/ M( _β^t ∂_β v^s+ _3^t ∂_3 v^s)
    - ∂_α w_α+
    1/M_l^t_l+.


The field components _α^b and w_α need to be eliminated.
From equation (E.16) we obtain


    _α^b=-1/ıω(  _αβ∂_β v^s+  _α 3∂_3  v^s+ıω/ M_α p^f -_α l_l).


Pre-multiplying all terms in  equation (E.15) by _α^t gives


    w_α =_α^t(ıω/η^f v^s -  1/η(_β∂_βp^f+_3∂_3p^f)+1/η_j f_j^f).


Using this in equation (<ref>) gives



    ∂_3  p^f    =   ıω_3^t^f v^s -η_3^t^-1_3 w_3
    
    -_3^t^-1_α_α^t(ıω^f v^s -  (_β∂_βp^f+_3∂_3p^f)+_j f_j^f)+ f_3^f,
          
    (1-_3^t^-1_α_α^t_3)∂_3  p^f    =   ıω_3^t^f v^s -η_3^t^-1_3 w_3
    
    -_3^t^-1_α_α^t(ıω^f v^s - _β∂_βp^f+_j f_j^f)+ f_3^f,
    ∂_3  p^f    =   1/(ıω_3^t^f v^s -η_3^t^-1_3 w_3
    -_3^t^-1_α_α^t(ıω^f v^s - _β∂_βp^f+_j f_j^f)+ f_3^f)
       =   1/(_3^t^-1_α_α^t_β∂_βp^f
    +ıω_3^t( I-^-1_α_α^t)^f v^s-η_3^t^-1_3 w_3
    
    -_3^t^-1_α_α^t_jf_j^f+f_3^f),



with


    =̱1-_3^t^-1_α_α^t_3.


Substituting  this into equation (<ref>) gives



    w_α   =   _α^t(ıω/η^f v^s -  1/η_β∂_βp^f+1/η_j f_j^f)
          
    -1/η_α^t_3(_3^t^-1_γ_γ^t_β∂_βp^f
    +ıω_3^t( I-^-1_γ_γ^t)^f v^s-η_3^t^-1_3 w_3
    
    -_3^t^-1_γ_γ^t_jf_j^f+f_3^f),
    
    w_α   =   -1/η(_α^t_β + 1/_α^t_3_3^t^-1_γ_γ^t_β)∂_β p^f
    +ıω/η_α^t( I-1/_3_3^t( I-^-1_γ_γ^t))^f v^s
          +1/_α^t_3_3^t^-1_3 w_3+1/η_α^t( I+
    1/_3_3^t^-1_γ_γ^t)_jf_j^f-1/η_α^t_3f_3^f.



We are now ready to eliminate _α^b and w_α from equations (<ref>) - (<ref>).
The expression for ∂_3 v^s, equation (<ref>), already has the desired form.
Substituting equations (<ref>) and (<ref>) into equation (<ref>), we obtain



    -∂_3_3^b    =   
    -ıω/η^f_α(_α^t_β + 1/_α^t_3_3^t^-1_γ_γ^t_β)∂_β p^f
    
    +ıω^b v^s-ω^2/η^f_α_α^t( I-1/_3_3^t( I-^-1_γ_γ^t))^f v^s
    
          +ıω^f( I+1/_α_α^t_3_3^t^-1)_3 w_3
    +ıω/η^f_α_α^t( I+1/_3_3^t^-1_γ_γ^t)_jf_j^f
    -ıω/η^f_α_α^t_3f_3^f
    
          -1/ıω∂_α(
     _αβ∂_β v^s + _α 3∂_3  v^s
    +ıω/ M_α p^f) +  f^b
    +1/ıω∂_α(
    _α l_l).



Upon substitution of equation (<ref>) and using


    ^-1_γ_γ^t   =   ^-1( I-_3_3^t)=  I-^-1_3_3^t,
    _3^t^-1_γ_γ^t_β   =   
    
    _3^t( I-^-1_3_3^t)_β=
     -_3^t^-1_3_3^t_β,


we obtain



    -∂_3_3^b    =   ∂_α( _α 3_33^-1_3^b)
    -ıω/η^f_α(_α^t_β - 1/_α^t_3_3^t^-1_3_3^t_β)∂_β p^f
    -1/ıω∂_α(ıω/ M_α p^f + _αβ∂_β v^s)
          +ıω^b v^s-ω^2/η^f_α_α^t( I-1/_3_3^t^-1_3_3^t)^f v^s
    
    +ıω^f( I+1/_α_α^t_3_3^t^-1)_3 w_3
    
          
    +ıω/η^f_α_α^t( I+1/_3_3^t^-1_γ_γ^t)_jf_j^f
    -ıω/η^f_α_α^t_3f_3^f
    
    +  f^b
    +1/ıω∂_α(
    _α l_l),



with 


    _α l   =   _α l-
     _α 3_33^-1_3l,
    _l   =   _l-  _l3_33^-1_3,


where _α 3= O,  _3= 0 and


    _αβ^t= _βα


on account of _jl^t=_lj.
Using equations (<ref>) and  (<ref>)  in equation (<ref>), we obtain


    ∂_3 p^f
       =   1/(-_3^t^-1_3_3^t_β∂_βp^f
    +ıω_3^t^-1_3_3^t^f v^s-η_3^t^-1_3 w_3
    
    -_3^t^-1_α_α^t_jf_j^f+f_3^f),


with 



    =̱1-_3^t^-1_α_α^t_3=
    _3^t( I-^-1_α_α^t)_3=
    _3^t^-1_3_3^t_3=
    (_3^t^-1_3)(_3^t_3)=
    (_3^t_3)(_3^t^-1_3)=
    _3^t_3_3^t^-1_3.


Substituting equations (<ref>) and (<ref>)
into equation (<ref>), we obtain


    ∂_3  w_3    =   ıω/ M p^f
      -  1/ M_β^t ∂_β v^s  
    
      + 1/ M_3^t  _33^-1(ıω_3^b
    +ıω/ M_3 p^f + _3β∂_β v^s-_3l_l
    )
           -∂_α(-1/η(_α^t_β - 1/_α^t_3_3^t^-1_3_3^t_β)∂_β p^f
    +ıω/η_α^t( I-1/_3_3^t^-1_3_3^t)^f v^s
          +1/_α^t_3_3^t^-1_3 w_3+1/η_α^t( I+
    1/_3_3^t^-1_γ_γ^t)_jf_j^f-1/η_α^t_3f_3^f )+
    1/M_l^t_l+,


or


    ∂_3  w_3    =   ıω/ M_3^t  _33^-1_3^b+ ıω/ M p^f + ıω/M^2_3^t  _33^-1_3 p^f
    +∂_α1/η(_α^t_β - 1/_α^t_3_3^t^-1_3_3^t_β)∂_β p^f
            -  1/ M_β^t ∂_β v^s
    -∂_αıω/η_α^t( I-1/_3_3^t^-1_3_3^t)^f v^s
          -∂_α(1/_α^t_3_3^t^-1_3 w_3+1/η_α^t( I+
    1/_3_3^t^-1_γ_γ^t)_jf_j^f-1/η_α^t_3f_3^f )+
    1/M_α^t_α+.


Equations (<ref>), (<ref>), (<ref>) and (<ref>) can be cast  in the form of matrix-vector wave equation  (1) in the main paper,
with the  wave vector = ( x,ω) and
 the source vector = ( x,ω) defined as


    = [ - _3^b;    p^f;    v^s;    w_3;        ],    = [ ıω/η^f_α_α^t( I+1/_3_3^t^-1_γ_γ^t)_jf_j^f
    -ıω/η^f_α_α^t_3f_3^f
    +  f^b+1/ıω∂_α(
    _αβ_β);                                                       1/(-_3^t^-1_α_α^t_jf_j^f+f_3^f);                                                                           _33^-1_3l_l;               -∂_α(1/η_α^t( I+
    1/_3_3^t^-1_γ_γ^t)_jf_j^f-1/η_α^t_3f_3^f )+
    1/M_α^t_α+ ]


and the operator matrix =( x,ω) having the form of equation  (2) in the main paper, with


    _11   =   [ _11^11 _11^12; _11^21 _11^22 ],  _12 = [ _12^11 _12^12; _12^21 _12^22 ],
    _21   =   [ _21^11 _21^12; _21^21 _21^22 ],  _22 = [ _22^11 _22^12; _22^21 _22^22 ],


where, using equation (<ref>),


    _11^11   =    -∂_α_α 3_33^-1,
    _11^12   =   -ıω/η^f(_α_α^t - _α_α^t_3_3^t^-1_3_3^t)_β∂_β 
    -∂_α1/M_α
       =   
    -ıω/η^f(I̱-_3_3^t - _α_α^t_3_3^t^-1_3_3^t)_β∂_β 
    -∂_α1/M_α
       =   -ıω/η^f(I̱-(_3_3^t + _α_α^t)_3_3^t^-1_3_3^t)_β∂_β 
    -∂_α1/M_α
       =   -ıω/η^f( I-1/_3_3^t^-1_3_3^t)_β∂_β 
    -∂_α1/M_α,
    
    
    
    _11^21   =    0^t,
    _11^22   =   -1/_3^t^-1_3_3^t_β∂_β,
    
    
    
    _12^11   =   ıω^b-1/ıω∂_α_αβ∂_β
    -ω^2/η^f_α_α^t(I̱-_3_3^t^-1_3_3^t)^f
       =   ıω^b-1/ıω∂_α_αβ∂_β
    -ω^2/η^f_α_α^t(_β_β^t+_3_3^t-_3_3^t^-1_3_3^t)^f
       =   ıω^b-1/ıω∂_α_αβ∂_β
    -ω^2/η^f_α_α^t(_β_β^t-_3_3^t^-1_3_3^t( I-_3_3^t))^f
       =   ıω^b-1/ıω∂_α_αβ∂_β
    -ω^2/η^f_α_α^t(-1/_3_3^t^-1_3_3^t)_β_β^t^f,
    _12^12   =   ıω/^f(_3+̱_α_α^t_3_3^t^-1_3 )
    =ıω/^f(_3_3^t+_α_α^t)_3_3^t^-1_3
       =   ıω/^f_3_3^t^-1_3,
    _12^21   =   ıω/_3^t^-1_3_3^t^f,
    _12^22   =   -η/_3^t^-1_3,
    
    
    
    _21^11   =   ıω_33^-1,
    _21^12   =   -ıω/M_33^-1_3,
    _21^21   =   -ıω/M_3^t_33^-1,
    _21^22   =   ıω/M^2_3^t_33^-1_3 +ıω/M
    +∂_α1/η(_α^t_β - 1/_α^t_3_3^t^-1_3_3^t_β)∂_β,



    _22^11   =   -_33^-1_3β∂_β,
    _22^12   =    0,
    _22^21   =    -  1/ M_β^t ∂_β
    -∂_αıω/η_α^t( I-1/_3_3^t^-1_3_3^t)^f,
    _22^22   =   -∂_α1/_α^t_3_3^t^-1_3.


For the special case of an isotropic medium we have


    _̧ijkl   =   ( K_ G-2/3 G_ fr)δ_ijδ_kl+ G_ fr(δ_ikδ_jl+δ_ilδ_jk),
    
     C_ij   =    Cδ_ij,
    ρ_ij^b   =   ρ^bδ_ij,
    ρ_ij^f   =   ρ^fδ_ij,
    
    k_ij   =   kδ_ij.


Here G_ fr
is the complex frequency-dependent shear modulus of the framework of the grains when the fluid is absent. 
The elastic parameters K_ G (Gassmann modulus), C and M are given by <cit.> 


    K_ G   =    K_ fr+ϕ  K^f+(1+ϕ) K^s Δ/1+ Δ,
    
     C    =    K^f+ K^s Δ/1+ Δ,
    
    
    
     M    =   1/ϕ K^f/1+ Δ,
    Δ   =    K^f/ϕ ( K^s)^2((1-ϕ)
    K^s- K_ fr),


where K^s and K^f are the solid and fluid  compression moduli
and K_ fr is the
compression modulus of the framework of the grains. These parameters can be
expressed in terms of Biot's parameters A, N, Q and R <cit.>, according to


    K_ G-2/3 G_ fr= A+2 Q+ R,  
     G_ fr= N,     C= Q+ R/ϕ,
        M= R/ϕ^2.


Hence,


    (_jl)_ik   =   _ijkl=_̧ijkl-C^2/Mδ_ijδ_kl=Sδ_ijδ_kl+N(δ_ikδ_jl+δ_ilδ_jk),
    
    S   =   A-Q^2/R,
    _j   =   C_j,
    ^b   =   ρ^b I,
    ^f   =   ρ^f I,
       =   k I,
    =   1.


With these substitutions, we obtain for the  source vector


    = [ ıωρ^f k/η_α f_α^f
    +  f^b+1/ıω∂_α(_αβ_β);                                   f_3^f;                             _33^-1_3l_l;              -∂_α(k/ηf_α^f)+
    1/M_α^t_α+ ]


and for the operator matrices


    _11^11   =    -∂_α_α 3_33^-1,
    _11^12   =   -ıωρ^f k/η_β∂_β-
    ∂_α1/M_α,
    _11^21   =    0^t,
    _11^22   =   0,
    
    
    
    _12^11   =    -1/ıω∂_α_αβ∂_β+ıω(ρ^b I
    +ıω(ρ^f)^2k/η_α_α^t),
    _12^12   =   ıωρ^f_3,
    _12^21   =   ıωρ^f_3^t,
    _12^22   =   -η/ k,


    _21^11   =   ıω_33^-1,
    _21^12   =   -ıωC/M_33^-1_3,
    _21^21   =   -ıωC/M_3^t_33^-1,
    _21^22   =   ıωC^2/M^2_3^t_33^-1_3 +ıω/M+∂_α k/η∂_α,
    
    
    
    _22^11   =   -_33^-1_3β∂_β,
    _22^12   =    0,
    _22^21   =    -ıω∂_αρ^f  k/η_α^t-1/M_β^t∂_β,
    _22^22   =   0.


Using this in equations (<ref>) and (<ref>) we obtain


    _11   =   [                       0                       0              -∂_1 S/K_c  ρ^f/ρ^E∂_1-∂_12CN/MK_c;                       0                       0              -∂_2
    S/K_c ρ^f/ρ^E∂_2-∂_2
    2CN/MK_c;                    -∂_1                    -∂_2                       0                       0;                       0                       0                       0                       0 ],
    _12   =   [ ıω(ρ^b-(ρ^f)^2/ρ^E)
    -1/ıω( ∂_1ν_1∂_1+
    ∂_2N∂_2)                  -1/ıω(∂_2 N ∂_1+ ∂_1 ν_2∂_2 )                                              0                                              0;                   -1/ıω(∂_2 ν_2∂_1+ ∂_1 N∂_2 )  ıω(ρ^b-(ρ^f)^2/ρ^E)-1/ıω(
    ∂_1N∂_1+ ∂_2ν_1∂_2)                                              0                                              0;                                              0                                              0                                          ıωρ^b                                          ıωρ^f;                                              0                                              0                                          ıωρ^f                                          ıωρ^E ],
          
    _21   =   [                              ıω/N                                 0                                 0                                 0;                                 0                              ıω/N                                 0                                 0;                                 0                                 0                            ıω/K_c                        -ıω C/MK_c;                                 0                                 0                        -ıω C/MK_c ıω(C^2/M^2K_c +1/M)-∂_α1/ıωρ^E∂_α ],
    _22   =   [                       0                       0                    -∂_1                       0;                       0                       0                    -∂_2                       0;               -S/K_c∂_1               -S/K_c∂_2                       0                       0; ∂_1ρ^f/ρ^E
    -2CN/MK_c∂_1 ∂_2ρ^f/ρ^E
    -2CN/MK_c∂_2                       0                       0 ],


where


    ρ^E   =   ρ^E( x,ω)=-η/ıω k,
    ν_1   =   ν_1( x,ω)=
    4N(S+N/K_c),
    ν_2   =   ν_2( x,ω)=
    2N(S/K_c),
    
     K_c   =    K_c( x,ω)= S+2 N.






§ PIEZOELECTRIC WAVES




Equations (F.9) - (F.14) in Appendix F in the main paper form the starting point for deriving a piezoelectric
matrix-vector wave equation 
for the quantities _3, E_0, v and H_0.
We separate the derivatives in the x_3-direction from the lateral derivatives, according to


    -∂_3_3    =   ıω v + ∂_α_α+  f,
    ∂_3  E_0    =   ıω_1  H_0 +ıω_3H_3+ E_3-   J_0^m,
    ∂_3  v   =   _33^-1(-ıω_3 -_3β∂_β v -ıω_3l( D_1l E_0+ D_3lE_3)+_3l h_l),
    ∂_3  H_0   =   ıω_1   E_0+ıω_3 E_3+  H_3 + ıω D_1k^t_k  -  J_0^e.


The field components  _α,  E_3 and H_3  need to be eliminated. We start by deriving an explicit expression for E_3. 
From equations (F.10) and (F.14) we obtain


    E_3   =    E_33^-1(- _3^t  E_0  +1/ıω^t  H_0
    
    + D_3k^t(1/ıω_k l∂_l  v  +_k l( D_1l E_0+ D_3lE_3)- 1/ıω_k l h_l)+1/ıω J_3^e),
          
    
    ( E_33- D_3k^t_k l D_3l)E_3   =   - _3^t  E_0  +1/ıω^t  H_0
    
    + D_3k^t(1/ıω_k l∂_l  v  +_k l D_1l E_0- 1/ıω_k l h_l)+1/ıω J_3^e,
    
    E_3   =   ()^-1(- ()^t  E_0  +1/ıω^t  H_0
    
    +1/ıω D_3k^t(_k β∂_β v +_k 3∂_3  v - _k l h_l)+1/ıω J_3^e),


with


    =    E_33- D_3k^t_k l D_3l,
    
    ()^t   =   _3^t - D_3k^t_k l D_1l.


The term ∂_3  v is eliminated from equation (<ref>) by substituting
equation (<ref>), hence 


    E_3   =   ()^-1(- ()^t  E_0  +1/ıω^t  H_0 +1/ıω D_3k^t(_k β∂_β v  - _k l h_l)+1/ıω J_3^e
          +1/ıω D_3k^t_k3_33^-1(-ıω_3 -_3β∂_β v -ıω_3l( D_1l E_0+ D_3lE_3)+_3l h_l)),
    
    (+ D_3k^t_k3_33^-1_3l D_3l)E_3   =   
    - ()^t  E_0  +1/ıω^t  H_0 +1/ıω D_3k^t(_k β∂_β v  - _k l h_l)+1/ıω J_3^e
          +1/ıω D_3k^t_k3_33^-1(-ıω_3 -_3β∂_β v -ıω_3l D_1l E_0+_3l h_l),
    
    E_3   =   ()^-1(- ()^t  E_0  +1/ıω^t  H_0+1/ıω D_3k^t(_k β∂_β v - _k l h_l)
    
    - D_3k^t_k 3_33^-1_3+1/ıω J_3^e),


with


    =   + D_3k^t_k 3_33^-1_3l D_3l= E_33- D_3k^t_kl D_3l,
    
    ()^t   =   ()^t+ D_3k^t_k 3_33^-1_3l D_1l=_3^t- D_3k^t_kl D_1l,
    _kl   =   _kl-_k 3_33^-1_3l,


where _k 3=_3 l= O and


    _αβ^t= _βα


on account of _jl^t=_lj.
Next, we derive an  expression for H_3 from equation (F.12), according to


    H_3=μ_33^-1(- _3^t  H_0  +1/ıω^t  E_0+1/ıω J_3^m).


From equation (F.14) we obtain the following expression for _α 


    -_α=1/ıω_αβ∂_β v+1/ıω_α 3∂_3  v+_α l( D_1l E_0+ D_3lE_3)-1/ıω_α l h_l,


from which ∂_3  v and E_3 need to be eliminated.
Substituting equation (<ref>) into equation (<ref>) yields


    ∂_3  v   =   _33^-1(-ıω_3 -_3β∂_β v -ıω_3l D_1l E_0+_3l h_l
          -ıω()^-1_3l D_3l(- ()^t  E_0  +1/ıω^t  H_0+1/ıω D_3k^t(_k β∂_β v - _k m h_m)
    
    - D_3k^t_k 3_33^-1_3+1/ıω J_3^e)),
    ∂_3  v   =   _33^-1(-ıω I'_3
    -ıω_3l D_1l' E_0 -_3β'∂_β v 
    -()^-1_3l D_3l^t  H_0
    
    -()^-1_3l D_3l J_3^e+_3m' h_m),


with


    I'   =    I-()^-1_3l D_3l D_3k^t_k 3_33^-1
    _3m'   =   _3m+ ()^-1_3l D_3l D_3k^t_k m,     m=β(=1,2), 3,
     D_1l'    =    D_1l  - ()^-1 D_3l ()^t.


Substituting equations (<ref>) and (<ref>) into equation (<ref>) gives


    -_α   =   1/ıω_αβ∂_β v
    +1/ıω_α 3_33^-1(-ıω I'_3
    -ıω_3l D_1l' E_0 -_3β'∂_β v 
    -()^-1_3l D_3l^t  H_0
          -()^-1_3l D_3l J_3^e+_3m' h_m) +_α l D_1l E_0
          +()^-1_α l D_3l(- ()^t  E_0  +1/ıω^t  H_0+1/ıω D_3k^t(_k β∂_β v - _k m h_m)
    
    - D_3k^t_k 3_33^-1_3+1/ıω J_3^e)- 1/ıω_α m h_m,
    
    -_α   =    - _α_3 - _α E_0 + 1/ıω_αβ' ∂_β v  - 1/ıω_α^t  H_0
    -1/ıω_α m' h_m -1/ıω_α J_3^e,


with


    _α   =   _α 3_33^-1 I'+()^-1_α l D_3l D_3k^t_k 3_33^-1=(_α 3+()^-1_α l D_3l D_3k^t_k 3)_33^-1,
    _α   =   _α 3_33^-1_3l D_1l'-_α l D_1l+()^-1_α l D_3l ()^t=
    -_α l( D_1l-()^-1 D_3l ()^t), 
    _α m'   =   _α m-_α 3_33^-1_3m'+()^-1_α l D_3l D_3k^t_km =_α m+()^-1_α l D_3l D_3k^t _k m,
    _α   =   ()^-1_α 3_33^-1_3l D_3l-()^-1_α l D_3l=-()^-1_α l D_3l,


where _α 3'= O and


    (_αβ')^t=_βα'.


We are now ready to eliminate  _α,  E_3 and H_3 from equations (<ref>) - (<ref>).
The expression for ∂_3 v, equation (<ref>), already has the desired form.
Substituting equation (<ref>) into (<ref>) gives


    -∂_3_3    =   ıω v + ∂_α(_α_3 + _α E_0 - 1/ıω_αβ' ∂_β v  + 1/ıω_α^t  H_0
    +1/ıω_αβ' h_β +1/ıω_α J_3^e)+  f
       =   ∂_α_α_3+∂_α_α E_0 +(ıω-1/ıω∂_α_αβ' ∂_β) v 
    +1/ıω∂_α_α^t  H_0 +1/ıω∂_α_αβ' h_β  +1/ıω∂_α_α J_3^e+ f.


Substituting equations  (<ref>) and (<ref>) into (<ref>) yields 


    ∂_3  E_0    =   ıω_1  H_0 +ıω_3μ_33^-1(- _3^t  H_0  +1/ıω^t  E_0+1/ıω J_3^m)-   J_0^m
          +(()^-1(- ()^t  E_0  +1/ıω^t  H_0+1/ıω D_3α^t(_αβ∂_β v - _αβ h_β)
    
    - D_3k^t_k 3_33^-1_3+1/ıω J_3^e))
       =   -()^-1 D_3k^t_k 3_33^-1_3+(_3μ_33^-1^t - ()^-1()^t ) E_0
    +1/ıω()^-1 D_3α^t_αβ∂_β v
          +(ıω(_1- _3μ_33^-1_3^t)+1/ıω()^-1^t)  H_0
    
    + _3μ_33^-1 J_3^m-   J_0^m-1/ıω()^-1 D_3α^t_αβ h_β+1/ıω()^-1 J_3^e.


Substituting equations  (F.14), (<ref>), (<ref>) and (<ref>)   into equation (<ref>) gives


    ∂_3  H_0   =   ıω_1   E_0+
     ıω_3 E_3
    
     +(μ_33^-1(- _3^t  H_0  +1/ıω^t  E_0+1/ıω J_3^m))
     
           
     - D_1k^t(_kβ∂_β v+_k3∂_3  v+ıω_kl( D_1l E_0+ D_3lE_3)- _kl h_l) 
      -  J_0^e
       =   ıω_1   E_0 + ıω()^-1(_3- D_1m^t_ml D_3l) (- ()^t  E_0  +1/ıω^t  H_0+1/ıω D_3k^t(_k β∂_β v - _k l h_l)
    
          - D_3k^t_k 3_33^-1_3+1/ıω J_3^e)
     +(μ_33^-1(- _3^t  H_0  +1/ıω^t  E_0+1/ıω J_3^m))
            - D_1k^t(_kβ∂_β v+_k3_33^-1(-ıω I'_3
    -ıω_3l D_1l' E_0 -_3β'∂_β v 
    -()^-1_3l D_3l^t  H_0
          -()^-1_3l D_3l J_3^e+_3m' h_m  )+ıω_kl D_1l E_0- _kl h_l) 
      -  J_0^e
       =   ıω( D_1k^t_k3_33^-1 I' -()^-1 D_3k^t_k 3_33^-1) _3  
     +(ıω( - ()^-1 ()^t )+1/ıωμ_33^-1^t)   E_0 
          +(()^-1 D_3k^t_k β  -  D_1k^t(_k β-_k3_33^-1_3 β') )∂_β v 
     +(()^-1^t -μ_33^-1_3^t )  H_0   
          -(()^-1 D_3k^t_km  -  D_1k^t(_km-_k3_33^-1_3m') ) h_m 
     + ()^-1 J_3^e+1/ıωμ_33^-1J_3^m-  J_0^e,


with


    =   _1 - D_1k^t(_k l D_1l-_k3_33^-1_3 l D_1l')=_1 - D_1k^t_k l D_1l- ()^-1 D_1k^t_k3_33^-1_3 l D_3l ()^t.


Using the following relations 

    D_1k^t_k3_33^-1 I' -()^-1 D_3k^t_k 3_33^-1   =   ( D_1k^t-()^-1 D_3k^t)_k 3_33^-1,
     - ()^-1 ()^t   =   _1 - D_1k^t_k l D_1l-()^-1()^t,
    
    ()^-1 D_3k^t_k m  -  D_1k^t(_k m-_k3_33^-1_3 m')    =   -( D_1k^t-()^-1 D_3k^t)_k m,


equation (<ref>) can be rewritten as


    ∂_3  H_0   =   ıω( D_1k^t-()^-1 D_3k^t)_k 3_33^-1_3
      + (ıω(_1 - D_1α^t_αβ D_1β-()^-1()^t)+1/ıωμ_33^-1^t)  E_0
          -( D_1α^t-()^-1 D_3α^t)_αβ∂_β v 
      +(()^-1^t -μ_33^-1_3^t )  H_0
          +( D_1α^t-()^-1 D_3α^t)_αβ h_β   
    + ()^-1 J_3^e+1/ıωμ_33^-1J_3^m-  J_0^e.


Equations  (<ref>), (<ref>), (<ref>) and (<ref>) can be cast in the form of matrix-vector wave equation  (1) in the main paper,

with the wave vector = ( x,ω) and 
 the source vector = ( x,ω) defined as


    = [ - _3;  E_0;    v;  H_0 ],    = [            f+1/ıω∂_α( _αβ' h_β)  - 1/ıω∂_α(()^-1_αβ D_3β J_3^e);      1/ıω(()^-1(J_3^e- D_3α^t_αβ h_β))-   J_0^m+_3μ_33^-1 J_3^m;                            _33^-1(_3m' h_m-()^-1_3l D_3l J_3^e); -  J_0^e   
    + ()^-1 J_3^e+1/ıω(μ_33^-1J_3^m)+ ( D_1α')^t_αβ h_β ]


and the  operator matrix =( x,ω) having the form defined in equation  (2) in the main paper, with


    _11   =   [ _11^11 _11^12; _11^21 _11^22 ],  _12 = [ _12^11 _12^12; _12^21 _12^22 ],
    _21   =   [ _21^11 _21^12; _21^21 _21^22 ],  _22 = [ _22^11 _22^12; _22^21 _22^22 ],


where


    _11^11   =   -∂_α_α=
    -∂_α(_α 3+()^-1_αβ D_3β D_3k^t_k 3)_33^-1,
    _11^12   =   ∂_α_α=-∂_α_αβ( D_1β-()^-1 D_3β ()^t),
    _11^21   =   ()^-1 D_3k^t_k 3_33^-1,
    _11^22   =   _3μ_33^-1^t - ()^-1()^t ,
    
    
    
    _12^11   =   ıω-1/ıω∂_α_αβ' ∂_β, 
    _12^12   =   1/ıω∂_α_α^t=-1/ıω∂_α()^-1_αβ D_3β^t,
    _12^21   =   1/ıω()^-1 D_3α^t_αβ∂_β,
    _12^22   =   ıω(_1- _3μ_33^-1_3^t)+1/ıω()^-1^t,
    
    
    
    _21^11   =   ıω(_33^-1 - ()^-1_33^-1_3l D_3l D_3k^t_k 3_33^-1),
    _21^12   =   -ıω_33^-1_3l D_1l' =-ıω_33^-1_3l( D_1l  - ()^-1 D_3l ()^t),
    _21^21   =   -ıω( D_1k')^t_k 3_33^-1=-ıω( D_1k^t-()^-1 D_3k^t)_k 3_33^-1,
    _21^22   =   ıω(_1 -()^-1()^t- D_1α^t_αβ D_1β)+1/ıωμ_33^-1^t,
    
    
    
    _22^11   =    - _33^-1_3β'∂_β=- _33^-1(_3β+ ()^-1_3l D_3l D_3α^t_αβ)∂_β, 
    _22^12   =    - ()^-1_33^-1_3l D_3l^t,
    _22^21   =   -( D_1α^t-()^-1 D_3α^t)_αβ∂_β,
    _22^22   =   ()^-1^t -μ_33^-1_3^t.







§ SEISMOELECTRIC WAVES



Equations (G.18) - (G.25) in Appendix G in the main paper form the starting point for deriving a seismoelectric
matrix-vector wave equation 
for the quantities _3^b, p^f, E_0,  v^s, w_3 and H_0 in vector . 
We separate the derivatives in the x_3-direction from the lateral
derivatives, according to


    -∂_3_3^b    =   ıωρ^b v^s
    +ıωρ^f(_α w_α+_3 w_3)
    +∂_α_α^b+
     f^b,
    
    
    
    
    ∂_3  p^f    =    ıωρ^f_3^t v^s-η/k(w_3-L E_3)+ f_3^f,
    
    
    
    ∂_3  E_0    =   ıωμ H_0 +
     E_3-
      J_0^m,
    
    
    
    ∂_3  v^s   =   -  _33^-1(ıω_3^b
    +ıω C/ M_3 p^f + _3β∂_β v^s
    - _3l_l
    ),



    ∂_3  w_3    =   ıω/ Mp -  C/ M( _β^t ∂_β v^s+ _3^t ∂_3 v^s)
    -∂_β  w_β+C/M_l^t _l+,
    
    
    
    
     ∂_3  H_0   =   ıω E E_0 -η/kL_α w_α
    +  H_3-
     J_0^e,


with 2× 1 unit vector (_α)_β=δ_αβ. The field components  _α^b,  w_α, E_3 and H_3 need to be eliminated.
Using equations (G.19) and (G.20),
we can eliminate the terms _α w_α and
∂_α_α^b from
equation (<ref>), yielding


    -∂_3_3^b    =   ıωρ^b v^s
    -ıωρ^f_α( -ıωρ^fk/η_α^t v^s +k/η∂_αp^f -L_α^t E_0)
    +ıωρ^f_3 w_3
    
          -1/ıω∂_α(
     _αβ∂_β v^s + _α 3∂_3  v^s
    +ıω C/ M_α p^f) +  f^b+ıωρ^f k/η_α 
    f_α^f+1/ıω∂_α(_α l_l),


or, upon substitution of equation (<ref>),


    -∂_3_3^b    =   ∂_α( _α 3_33^-1_3^b)
    -ıωρ^f k/η_α∂_αp^f
    -1/ıω∂_α(
    ıω/ M_α p^f + _αβ∂_β v^s
    )
          +ıω( ρ^b I_3
    +ıω( ρ^f)^2k/η_α_α^t)  v^s
    +ıωρ^f_3 w_3 +ıωρ^fL_α_α^t E_0 
    
    +  f^b+ıωρ^fk/η_α f_α^f+1/ıω∂_α(_α l_l),


with  I_3 being a 3× 3 identity matrix and


    _α l   =   _α l- _α 3_33^-1_3l,
    _l   =   C(_l-  _l3_33^-1_3),


where _α 3= O,  _3= 0 and


    _αβ^t= _βα


on account of _jl^t= _lj.
Using equation (G.23), we eliminate E_3 from equation
(<ref>), yielding


    ∂_3  p^f    =   ıωρ^f_3^t v^s-η/ k(1-1/ıω Eη/ kL^2)
     w_3
    
    +1/ıω Eη/ kL^t
      H_0
    +1/ıω Eη/ kL
    J_3^e+ f_3^f.


Using equation (G.23), we eliminate E_3 from
equation (<ref>), yielding


    ∂_3  E_0    =   (1/ıω Eη/ kL w_3)
    +ıωμ H_0 +
    (1/ıω E^t  H_0)
    
    -  J_0^m
    +(1/ıω E J_3^e).


Using equations (<ref>) and (G.19), we
eliminate  the terms ∂_3 v^s and ∂_βw_β from
equation (<ref>), according to


    ∂_3  w_3    =    C/ M_3^t _33^-1(ıω_3^b
    +ıω C/ M_3 p^f )+ ıω/ M p^f+∂_β( k/η∂_βp^f -L_β^t E_0
    -ıωρ^fk/η_β^t v^s )
            -  1/ M_β^t ∂_β v^s
    ++1/M_α^t_α
    -∂_β(  k/η f_β^f).



Using equations (G.19) and (G.25), we eliminate
 w_α and H_3 from equation (<ref>),
yielding


    ∂_3  H_0   =   L_α(  ∂_αp^f
    -η/ kL_α^t E_0 -ıωρ^f_α^t v^s )
    
    +ıω E E_0
    +(1/ıωμ^t
      E_0) -  J_0^e +(1/ıωμ J_3^m)-L_α f_α^f.


Equations (<ref>), (<ref>), (<ref>),  (<ref>),
 (<ref>) and (<ref>) can be cast in the form of matrix-vector wave equation  (1) in the main paper, 
 with the wave vector =(,ω) 
 and the source vector = ( x,ω)  defined as


    = [ - _3^b;    p^f;    E_0;    v^s;    w_3;    H_0 ],    = [ f^b+ıωρ^fk/η_α f_α^f+
    1/ıω∂_α(_αβ_β);            1/ıω Eη/ kL
    J_3^e
    + f_3^f;            -  J_0^m+
    (1/ıω E
     J_3^e);                          _33^-1_3l_l;          +1/M_α^t_α-∂_β(  k/η f_β^f);    -  J_0^e +(1/ıωμ J_3^m)-L_α
    f_α^f ]


and the operator matrix =( x,ω) having the form of equation  (2) in the main paper, with


    _11   =   [ _11^11 _11^12 _11^13; _11^21 _11^22 _11^23; _11^31 _11^32 _11^33 ],  _12 = [ _12^11 _12^12 _12^13; _12^21 _12^22 _12^23; _12^31 _12^32 _12^33 ],
    _21   =   [ _21^11 _21^12 _21^13; _21^21 _21^22 _21^23; _21^31 _21^32 _21^33 ],  _22 = [ _22^11 _22^12 _22^13; _22^21 _22^22 _22^23; _22^31 _22^32 _22^33 ],


where


    _11^11   =    -∂_α_α 3_33^-1,
    _11^12   =   -ıωρ^fk/η_α∂_α-
    ∂_α1/ M_α,
    _11^13   =   ıωρ^f
    L_α_α^t,
    
    
    
    _12^11   =   -1/ıω∂_α_αβ∂_β+ıω( ρ^b I_3
    +ıω( ρ^f)^2k/η_α_α^t),
    
    
    
    _12^12   =   ıωρ^f_3,
    _12^21   =   ıωρ^f_3^t,
    _12^22   =   -η/ k(1-1/ıω Eη/ kL^2),
    _12^23   =   1/ıω Eη/kL^t,
    _12^32   =   1/ıω Eη/ kL,
    _12^33   =   ıωμ I_2 +
    1/ıω E^t,
    
    
    
    _21^11   =   ıω_33^-1,
    _21^12   =   -ıω C/ M_33^-1_3,
    _21^21   =   -ıω C/ M_3^t _33^-1,
    _21^22   =   ıω C^2/ M^2_3^t _33^-1_3 +ıω/ M+∂_β k/η∂_β,
    
    
    
    _21^23   =   -∂_βL_β^t,
    _21^32   =   L_α∂_α,
    _21^33   =   ıω E I_2+
    1/ıωμ^t
    -η/ kL^2_α_α^t
    
    =(ıω E-η/ kL^2) I_2+
    1/ıωμ^t,
    
    
    
    _22^11   =   - _33^-1_3β∂_β,
    _22^21   =    -ıω∂_βρ^f k/η_β^t-1/ M_β^t
    ∂_β,
    _22^31   =   -ıωρ^fL_α_α^t,


where I_2 is a 2× 2 identity matrix. 
The submatrices that have not been listed here are zero.
Using this in equations (<ref>) and (<ref>), we obtain



    _11   =   [                           0                           0                -∂_1
     S/ K_c ρ^f/ρ^E∂_1-∂_1
    2 C N/ M K_c                      ıωρ^fL                           0;                           0                           0                -∂_2
     S/ K_c ρ^f/ρ^E∂_2-∂_2
    2 C N/ M K_c                           0                      ıωρ^fL;                        -∂_1                        -∂_2                           0                           0                           0                           0;                           0                           0                           0                           0                           0                           0;                           0                           0                           0                           0                           0                           0;                           0                           0                           0                           0                           0                           0 ],





    _12   =   (ıω( ρ^b-( ρ^f)^2/ρ^E)
    -1/ıω( ∂_1 ν_1∂_1+
    ∂_2 N∂_2)
        -1/ıω(∂_2  N∂_1+ ∂_1  ν_2∂_2 ) 
    
    
    -1/ıω(∂_2  ν_2∂_1+ ∂_1  N∂_2 )
       ıω( ρ^b-( ρ^f)^2/ρ^E)-1/ıω(
    ∂_1 N∂_1+ ∂_2 ν_1∂_2) 
    
     0    0 
    
     0     0 
    
     0     0  
    
     0     0 
    .
          
    
    
    
          .
    0     0   0   0
    
    0     0   0   0
    ıωρ^b    ıωρ^f   0   0 
    ıωρ^f    ıωρ^E(1+
     ρ^E/ EL^2)   ρ^E/ EL∂_1
       ρ^E/ EL∂_2
    
    0    -∂_1ρ^E/ EL   ıωμ-1/ıω∂_1
    1/ E∂_1    -1/ıω∂_1
    1/ E∂_2
    
    0     -∂_2ρ^E/ EL   -1/ıω∂_2
    1/ E∂_1
       ıωμ-1/ıω∂_2
    1/ E∂_2
    ),



    _21   =   [                                  ıω/ N                                      0                                      0                                      0                                      0                                      0;                                      0                                  ıω/ N                                      0                                      0                                      0                                      0;                                      0                                      0                                ıω/ K_c                          -ıω  C/ M K_c                                      0                                      0;                                      0                                      0                          -ıω  C/ M K_c ıω(  C^2/ M^2 K_c
    +1/ M)-∂_β1/ıωρ^E∂_β                                  -∂_1L                                  -∂_2L;                                      0                                      0                                      0                                   L∂_1                  ıωε-σ
    -1/ıω∂_2
    1/μ∂_2                         1/ıω∂_2
    1/μ∂_1;                                      0                                      0                                      0                                   L∂_2                         1/ıω∂_1
    1/μ∂_2                   ıωε-σ
    -1/ıω∂_11/μ∂_1 ],



    _22   =   [                           0                           0                        -∂_1                           0                           0                           0;                           0                           0                        -∂_2                           0                           0                           0;                 - S/ K_c∂_1                 - S/ K_c∂_2                           0                           0                           0                           0; ∂_1ρ^f/ρ^E
    -2 C N/ M K_c∂_1 ∂_2ρ^f/ρ^E
    -2 C N/ M K_c∂_2                           0                           0                           0                           0;                     -ıωρ^fL                           0                           0                           0                           0                           0;                           0                     -ıωρ^fL                           0                           0                           0                           0 ],


with ρ^E, ν_1, ν_2, S and K_c defined in section <ref>.




