Polarization effects of seismic waves on the basis of radiative transport theory

Summary

Radiative transfer theory provides a good framework for the study of multiple scattering in the randomly inhomogeneous lithosphere. Envelopes of high-frequency seismograms (mainly S coda waves) of local earthquakes have been synthesized on the basis of this theory, and inversions for some Earth parameters such as intrinsic attenuation, scattering attenuation and degree of non-isotropic scattering have been carried out. However, a scalar model has often been assumed because of its mathematical relative simplicity. The simplification amounts to neglecting the polarized nature of the underlying motion. This approach is only valid for long lapse times when S waves become unpolarized because of high-order scattering, and cannot be justified by only assuming that the source is unpolarized. We show that incoming unpolarized S waves can be up to 80 per cent polarized after single scattering. Depolarization of S waves after multiple scattering is studied by a Monte Carlo method. We show that the scattering of S waves off different kinds of inhomogeneities gives rise to different polarization and depolarization patterns. Consequently, polarization should provide valuable information for the understanding of the physics of wave motion and the properties of the Earth’s lithosphere.

coda, heterogeneity, Monte Carlo method, polarization, radiative transfer theory, scattering

1 Introduction

Coda wave analysis has attracted much interest in recent years. The coda is composed of scattered waves that sample the random medium in a similar way to the averaging of direct waves over a seismic network. Therefore, their study provides a good tool to understand the random inhomogeneities of the Earth’s lithosphere and estimate its properties.

Since the pioneering work of Aki (1969) on the seismic coda of local earthquakes, many authors (e.g. Aki & Chouet 1975; Sato 1977; Frankel & Wennerberg 1987; Herraiz & Espinosa 1987) have used phenomenological models to fit the total scattering coefficient and the coda wave attenuation to real S-coda data. More recently, theoretical (e.g. Wu 1985; Wu & Aki 1988; Shang & Gao 1988; Zeng et al. 1991; Sato 1993, 1994) and numerical studies (e.g. Gusev & Abubakirov 1987; Hoshiba 1991; Margerin et al. 1998) based on the radiative transfer theory (Chandrasekhar 1960; Ishimaru 1978) have been applied to synthesize coda envelopes. In all of those studies, the vector nature of the underlying wave motion is not taken into account, and a simplified scalar version of the motion is analysed. The simplification amounts to neglecting the polarized nature of shear waves. In addition, isotropic scattering is almost always assumed.

The aim of the present paper is to analyse the polarization effects of shear waves. Specifically, we follow the model derived in Ryzhik et al. (1996) for the propagation of elastic wave energy in random media to compute the spatio–temporal distributions of the energy density and of the degree of polarization. Starting from the full linear elastic wave equations, they derive a system of coupled radiative transfer equations for the angularly resolved P-wave energy density and S-wave coherence matrix. Based on the observation that S coda waves are composed primarily of S waves (Tsujiura 1978), we do not take into account the coupling between P and S waves. These observations are supported by the theory developed by Aki (1992), Zeng (1993) and Papanicolaou et al. (1996a,b), where the ratio of the energy conversion coefficients as functions of the P- to S-wave velocity ratio γ is given. In practice, γ≈√ 3. The scattered energy conversion from P to S waves is 2γ⁴≈ 18 times larger than the S- to P-wave conversion. It is also shown in Sato (1984) that the ratio of the radiated S-wave energy to theP-wave energy from a point shear dislocation is 1.5 γ⁵≈ 23.4. Both factors lead to a preponderance of S waves.

The spatio–temporal distributions of the energy density and the degree of polarization are computed by a Monte Carlo method. Unlike scalar equations, matrix-valued equations do not have a direct probabilistic representation suitable for Monte Carlo simulation. We follow the recent probabilistic representation of polarized radiative transfer equations given in Bal et al. (2000) to present a simple implementation for elastic radiative transfer; see also Bal & Moscoso (2000) for an application to electromagnetic radiative transfer.

This paper considers media with two different kinds of random fluctuations: fluctuations in the density and fluctuations in the Lamé parameters. We show that the scattering laws corresponding to these different fluctuations provide different energy coda envelopes. We present the spatio–temporal distribution of the polarization and depolarization of S waves. We show that these effects are very important for local earthquakes, and that they disappear only for large lapse times. In addition, they are very sensitive to the scattering law.

This paper is organized as follows. We describe in Section 2 our model for S-wave propagation. We show in Section 3 how to solve it by the Monte Carlo method. We study the single-scattering approximation assuming initially unpolarized waves in Section 4. In Section 5 we compute the effects of multiple scattering in the free space using the Monte Carlo method. In Section 6 we deal with half-space and slab problems, and show the degree of polarization at the free surface. The last section contains our conclusions. We have also included three Appendices, where we present some details of the probability scattering density, the single-scattering approximation and the implementation of the Monte Carlo method.

2 Radiative Transfer Equations for Shear Waves

Shear waves, like electromagnetic waves, cannot be completely described by energy density, but require a 2 × 2 non-negative Hermitian matrix,

which accounts for their polarized structure (Ryzhik et al. 1996; Chandrasekhar 1960). The parameter I describes the total energy density of the shear waves, whereas Q, U and V describe its polarization. Fully polarized waves satisfy the relation I²=Q²+U²+V², or equivalently detW=0. The polarization of a fully polarized transverse wave can be characterized by the end point of an arrow representing the wave motion. In this case, the arrow describes an ellipse with principal axes in the directions e₁ and e₂ and making an angle χ with a chosen orthogonal system of reference (e_l, e_r) transverse to the wave motion (see Fig.1a). The Stokes parameter Q is defined as the difference between the energy density of the waves linearly polarized along the direction e_l, I_l, and along the direction e_r, I_r. The same holds for U but with reference to two directions rotated by 45° with respect to (e_l, e_r). Finally, V determines the difference between the energy densities that are - or left-hand circularly polarized. Then, we can write

1

. (a) Polarization ellipse for a transverse wave. (b) Coordinate system for the polarization parameters of the incident and scattered waves in single scattering. The scattering plane is defined by the incident and scattered directions, and ∥ and perp stand for parallel and perpendicular directions to the scattering plane, respectively, in planes transverse to the incident and scattered waves.

where β is the angle (between 0 and π/2) whose tangent is the ratio of the principal axes of the ellipse. Consequently, for linearly polarized waves V is equal to 0. When the end point of the vector moves completely randomly, the wave is unpolarized and Q²+U²+V²=0. However, as we will see, the waves are generally in intermediate states of polarization, i.e. partially polarized. In this case we have I²≥Q²+U²+V². A scalar measure to quantify the amount of polarized energy density to the total amount of energy density is given by the degree of polarization,

which takes values between 0 (for unpolarized waves) and 1 (for fully polarized waves).

In this paper we consider a homogeneous constant-velocity background medium with weak random inhomogeneities on the scale of the wavelength. Inhomogeneities are variations of the elastic coefficients and mass density due to changes in material composition or the presence of microcracks. Upon interacting with these weak fluctuations, waves with wave vector k′ at point x are scattered into any direction {k^ with wave vector k, where k^=k/|k|. We do not consider P-to-S conversion and dissipation. With these assumptions, energy conservation is expressed by the transport equation (Ryzhik et al. 1996),

where v is the velocity of the homogeneous background and ∇_x is the transpose vector of (∂/∂x, ∂/∂y, ∂/∂z).

The local speed of propagation of shear waves is given by v(x)=√μ(x)/ρ (x). However, complete information about the spatial dependence of the elastic constant μ and the density ρ is not available for the whole Earth lithosphere, and a statistical treatment is necessary. For statistically homogeneous and statistically isotropic media, Ryzhik et al. (1996) derived the differential scattering tensor in terms of the power spectral densities R^_μμ, R^_ρρ and R^_μρ of the inhomogeneities. The power spectral density of the fluctuations of the density ρ (x) is the Fourier transform of the two-point correlation function

where ξ_ρ is the fractional fluctuation of the density and 〈˙〉 stands for ensemble averaging over realizations. Eq.(5) is independent of x because the fluctuations are assumed to be stationary (see Ryzhik et al. 1996). The power spectral densities R^_μμ and R^_μρ are given by analogue expressions. In this paper we assume R^_μρ=0, that is, there are no correlations between fluctuations in μ and ρ, The differential scattering tensor then has the form

where

are scalar functions due to fluctuations in the density ρ and the Lamé parameter μrespectively. Note that the delta functions in eq.(6) imply that the frequency ω=v|k| is unchanged by scattering. Therefore, different frequencies satisfy uncoupled equations. We consequently assume that the modulus |k| of the wavenumber of our initial condition W₀(x, k) is fixed.

The 2 × 2 matrices T(k, k′) and Γ(k, k′) are defined by

with the vectors (k^, z⁽¹⁾(k), z⁽²⁾(k)),

forming a local orthonormal basis in IR³. The differential scattering tensor (6) and the total scattering cross-section Σ(k) are related by

where II₂ is the 2 × 2 identity matrix, and

Here, cosΘ=k^˙k^′. Although expressions (10) and (11) can be checked by direct computation, we can argue the functional dependence as follows. The assumption of a statistically isotropic medium implies that the scattering coefficient depends only on k^˙k^′, leaving a trivial integration over the azimuthal angle. Moreover, assuming that R^_μμ=0 (respectively R^_ρρ=0), an unpolarized incident wave scattering in a direction making an angle Θ with the direction of incidence becomes partially plane polarized with a ratio of the energy densities equal to 1 /cos²Θ [respectively (cos²Θ−sin²Θ)²/cos²Θ] in directions perpendicular and parallel to the plane of scattering, which is spanned by the incident and scattered directions (see Fig.1b). This feature is important in analysing the single-scattering approximation.

Eq.(4) is a good approximation of the wave energy transport when the dominant wavelength is short compared to the macroscopic properties of the medium and the correlation length of the inhomogeneities is of the same order as the wavelength. In addition, fluctuations of the inhomogeneities must not be so strong as to cause localization. Let us emphasize that the coherence matrix W as well as the radiative transfer equation (4) are represented in a given orthonormal basis (e_x, e_y, e_z). The change of representation by rotation of the basis is by no means straightforward (Chandrasekhar 1960).

For the calculations that follow, it is convenient to render eq.(4) dimensionless by normalizing time, length and related quantities in terms of the S-wave velocity v and the mean free path 1/Σ [about 4kms^{− 1} and of the order of 100km, respectively, for frequencies of 1–30Hz (Sato & Fehler 1998)]. The initial condition that we consider in this paper is an unpolarized isotropic source term with frequency ωof the form W₀(x, k)=w₀II₂δ(x). We define the following adimensionalized quantities:

After integration over the modulus of the wave vector in eq.(4) using eq.(6), and dropping the tilde symbols for simplicity, we obtain the dimensionless equation

Here, as a simplification, we have assumed that R^_ρρ and R^_μμ are constant. Note that the real and imaginary parts of W are uncoupled since the matrices T and Γ are real. This means that the circularly polarized component of the waves evolves independently of the other components. Consequently, we assume without loss of generality that V=0, which corresponds to the case of linearly polarized shear waves. In the following we consider only the dimensionless radiative transfer equation (13).

To sum up, shear waves cannot be accurately described by a scalar quantity (i.e. the energy density). They require the intro-duction of the 2 × 2 coherence matrix (1), which accounts for their polarized structure. This paper analyses the polarization and depolarization effects of S waves. To simplify the analysis, we consider two examples: (i) fluctuating Lamé parameter μ and constant density ρ, and (ii) constant Lamé parameter μ and fluctuating density ρ, We will refer to them as the first and second cases, respectively.

3 Monte Carlo Method

We solve the radiative transfer equation (13) by the Monte Carlo method. This method enables us to deal with complicated 3-D geometries with reasonable storage requirements and is relatively easy to implement. For these reasons, it is very popular in radiation transport, heat transfer and nuclear reactor theories (Kalos & Whitlock 1986; Spanier & Gelbard 1969). The radiative transfer equations for scalar waves have the same structure as transport equations for particles. The central characteristic of Monte Carlo methods is to simulate the behaviour of one particle at a time. The average over a large number of realizations of particle paths yields an estimate of the particle probability density.

The simulation of the particle trajectories is threefold: (i) description of the particle trajectories between successive scattering events, (ii) estimation of the time of the next scatter-ing with the underlying medium, and (iii) at the time when the particle undergoes scattering, computation of the outgoing velocity of the particle. The last two steps are random processes.

As it is, this method does not hold for coherence matrices, whose entries do not correspond only to energies, but also to polarization parameters. Moreover, the probability distribution of the outgoing velocity after scattering does not only depend on the incident direction, as for scalar particle transport, but also on its state of polarization.

Nonetheless, it is possible to adapt these ideas for S waves by considering the polarization parameters as independent variables, that is, each particle will be characterized by its position X(t), velocity K(t) and polarization state P(t). We refer to Bal et al. (2000) for the details. Let us summarize here the main ideas.

We rewrite the matrix W as

where

Note that the energy density I is the trace of W(t, x, k). In addition to the position X(t) and wave vector (velocity) K(t) we now have a pair of polarization processes p₁(t) and p₂(t) that describe the polarization parameters Q and U (see eq.15). The position X(t) is continuous in space, while the velocity K(t)and polarization parameters P(t) jump at the scattering times.

The evolution (X(t), K(t), P(t)) of each particle follows these rules. Between scatterings, trajectories evolve along straight characteristics defined by

At time s, the probability that the next scattering takes place after time t follows the Poisson distribution

When a scattering event takes place, we have to define the outgoing velocity and polarization parameters. The probability density for a jump from wave vector k and polarization matrix P into wave vector k′ and polarization matrix P′ is

where the matrix S stands for T or Γ (see eq.6), C is a normalization constant and Tr stands for the trace of the matrix. The matrix of the scattered polarization parameters P_s(k,P,k′), which is symmetric and has null trace, is determined by (Bal et al. 2000)

For the reader’s convenience we have postponed to Appendix A the computation of the probability measure π and the scattered polarization parameters. We note that in these expressions the jump probability does not depend on the outgoing polarization parameters P′ and that the outgoing polarization parameters are functions of the incoming and outgoing velocities k and k′, and of the incoming polarization parameters P.

The history of each particle is constructed as follows. The history (X(t), K(t), P(t)) starts at (x₀, k₀, P⁰). The time of first scattering τ₁ is given by the probability density (17). The process follows (16) between 0 ≤t<τ₁. At τ₁ the process jumps to a new velocity K(τ₁) and a new polarization state P(τ₁) according to (18) and (19), respectively. Then, a new scattering time τ₂ is obtained according to (17). The process is continued in an obvious manner until a final time T_f.

Using the above construction of the jump process (X(t), K(t), P(t)), we approximate its probability density function ρ (t, x, k, p₁, p₂)by the classical Monte Carlo methods (see Spanier & Gelbard 1969), provided the initial distribution ρ (0, x, k, p₁, p₂) is known. The main property of the augmented process is that the linear moments with respect to the polarization parameters of the probability density function ρ (t, x, k, p₁, p₂) solve the radiative transfer equation (13). In other words, we have

We readily verify that the following initial condition models an isotropic unpolarized local pulse located at x=0:

Some details on the implementation of the Monte Carlo method are given in Appendix C. The validity of our code was checked by comparing the numerical results with the analytical single-scattering solution computed in Section 4 and with the diffusion solution for large lapse times. The agreement was perfect in both cases.

4 Single scattering

The application of single-scattering models (e.g. Aki 1969; Sato 1977) is based on the assumption of weak scattering, for which the mean free path (equal to 1 in our adimensional units) is much longer than the source–receiver distance. In addition to being a fairly good approximation, they also provide a good understanding of how polarization is generated by the interaction of unpolarized waves with random inhomogeneities. Therefore, single scattering may prove extremely useful in the estimation of the characteristics of the Earth’s lithosphere from high-frequency seismic data.

Here we consider inhomogeneities of the Earth’s lithosphere due to fluctuations of the Lamé parameter μ and fluctuations of the density ρ, successively. We will refer to them as the first and second cases, respectively, and they will be distinguished by the superscripts μ and ρ, In this section we analyse the angular dependence of the scattering patterns and calculate the spatio–temporal distribution of the energy density and of the degree of polarization for the extreme case of unpolarized sources using the relative simplicity of the spherical geometry (see Fig.2). The generalization to partially polarized sources is straightforward.

2

. Single-scattering geometry.

For this purpose, it is convenient to use the energy densities I_∥ and I⊥ in directions parallel and perpendicular, respectively, to the plane of scattering (which contains the incident and scattered directions; see Fig.1b). We denote the incident energy by

since we have assumed that V⁽ⁱ⁾=0. The scattered energy density in the direction Θ is proportional to M^{μ, ρ}J, where M^μ and M^ρ denote the phase matrices for the fluctuations in the Lamé parameter and for the fluctuations in the density, respectively. In this basis they are given by

Figs3(a) and (b) show the angular pattern of the total scattered energy density I (solid lines) for interactions between incident unpolarized S waves and inhomogeneities in μ and ρ, respectively. The energy densities I_∥ and I_perp are also shown, as dotted and dashed lines, respectively. In the two cases, the total energy density scattered in the forward and backward directions is two times larger than the energy scattered at right angles with respect to the incident wave. However, the scattered energy density attains its minimum at Θ=45° for the first case and at Θ=90° for the second. The difference between the scattering patterns is clear. Figs3(c), (d) and (e), (f) are the corresponding angular patterns when the incident waves are partially polarized. In Figs3(c) and (e), I⁽ⁱ⁾_∥=0.25 I⁽ⁱ⁾ and I⁽ⁱ⁾⊥ =0.75 I⁽ⁱ⁾, and in Figs3(d) and (f), I⁽ⁱ⁾_∥=0.75 I⁽ⁱ⁾ and I⁽ⁱ⁾⊥=0.25 I⁽ⁱ⁾, where I⁽ⁱ⁾ is the total incident energy density. These figures clearly show that the scattering angular patterns also depend on the incoming polarization state of the wave.

3

. Angular pattern (as a function of Θ) of the scattered energy density. The total scattered energy density I=I_∥+I_perp is represented with solid lines, whilst the scattered energy density in the parallel direction I_∥ and in the orthogonal direction I_perp are represented with dotted and dashed lines, respectively. (a), (c) and (d) correspond to the interaction of shear waves with inhomogeneities in μ. (b), (e) and (f) correspond to the interaction of shear waves with inhomogeneities in ρ, In (a) and (b), the incident waves are unpolarized. They are partially polarized in (c)–(f). In (c) and (e), I⁽ⁱ⁾_∥=0.25 I⁽ⁱ⁾ and I⁽ⁱ⁾⊥=0.75 I⁽ⁱ⁾, where I⁽ⁱ⁾ is the total incident energy density. In (d) and (f), I⁽ⁱ⁾_∥=0.75 I⁽ⁱ⁾ and I⁽ⁱ⁾⊥=0.25 I⁽ⁱ⁾. The scatterer is located at the origin of each figure.

Let us now compute the spatio–temporal distribution of the energy density and of the degree of polarization for unpolarized sources (I⁽ⁱ⁾_∥=I⁽ⁱ⁾⊥ =I⁽ⁱ⁾/2 and U⁽ⁱ⁾=0). Our point source spherically radiates unpolarized energy, which is received at time t, after single scattering, at a point located at a distance r from the source. Note that after single scattering, U remains equal to 0, and the energy densities of the scattered wave in directions parallel and perpendicular, respectively, to the plane of scattering are in the ratio [cos²Θ−sin²Θ]² to cos²Θ for the first case and cos²Θ to 1 for the second.

To compute the energy density and the degree of polarization arriving with angle ϕ at a distance r from the source at time t, we need to know the relation between the scattering angle Θ, the distance r₁≡|r₁| where the scattering took place and the direction ϕ (see Fig.2). Some algebra yields

The degree of polarization α and the energy density I arriving with angle ϕ are given in the single-scattering approximation by

where r₁ and Θ are given by (23). The common factor in I^{μ, ρ} comes from the geometrical spread of the energy density from the source to the point r. Its derivation is explained in Appendix B. Note that (24) and (26) correspond to the definition (3), since U=V=0 and Q=I_perp−I_∥. The degrees of polarization in (24) and (26) measure the difference of energy density in the directions perpendicular and parallel to the scattering plane normalized by the total energy density. Its value is 0 for unpolarized scattered waves and 1 for fully polarized scattered waves.

In Figs4(a) and 5(a), we show the spatial distribution of the polarization after one scattering for different angles ϕ for the first and second cases, respectively. It is clearly see that the polarization pattern enables us to distinguish between the two cases. In the first case, we observe two or three maxima, depending on the angle ϕ, where the degree of polarization equals 1 (fully polarized wave), with minima in between, where the degree of polarization is close to zero (unpolarized wave). In the second case, the polarization is smoother, with only one maximum in the interval [0, 1]. The average degree of polarization,

4

. (a) Spatial distributions of the degree of polarization α^μ(r, ϕ, 1)(eq.24) for ϕ=π/2 (solid line), ϕ=π/3 (dashed line), ϕ=π/6 (dash-dotted line) and ϕ=π/8 (dotted line). (b) Spatial distribution of the averaged polarization α^μ₂(r, 1) (eq.28).

5

. (a) Spatial distributions of the degree of polarization α^ρ(r, ϕ, 1)(eq.26) for ϕ=π/2 (solid line), ϕ=π/3 (dashed line), ϕ=π/6 (dash-dotted line) and ϕ=π/8 (dotted line). (b) Spatial distribution of the average polarization α^ρ₂(r, 1) (eq.28).

is plotted in Figs4(b) and 5(b). Observe that α^{μ, ρ}₂ exactly vanishes at the source location. This is because the scattering angle is π for all directions ϕ and because unpolarized waves scattered in the back direction remain unpolarized. Also note the three peaks in Fig.4(b), while only one is present in Fig.5(b)

Although (28) represents an average over angles ϕ, we prefer to use

which represents the percentage of fully polarized diffuse energy. The new ratio α₁ is easier to simulate numerically than α₂ when multiple scattering is included (see Section 5). The reason is that the directions of highest energy densities, for which the resolution is better in our Monte Carlo simulations, are more represented in the average α₁ than in α₂.

We plot in Fig.6 the spatio–temporal distributions of α^μ₁ (solid line) and α^ρ₁ (dashed line). In both cases the maximal value of the degree of polarization is around 0.8. This means that up to 80 per cent of the energy is polarized after single scattering. The difference between both cases is again clear. In the first case (solid line), a first peak is reached very quickly (see Fig.6b), while fluctuations in the density ρ (dashed line) do not generate polarized waves at this time. The second peak coincides in both cases. We also observe that the degree of polarization α^ρ₁ decreases monotonously to zero, while α^μ₁ has a third peak. At large times, the percentage of fully polarized energy differs by a factor of greater than 2 between the two cases. In Appendix B, we have derived analytical expressions for I^ρ(r, t), I^μ(r, t) and α^ρ₁(r, t).

6

. (a) Spatial distributions at t=1 and (b) time traces at r=1 of degree of polarization (29). Solid lines correspond to fluctuations in μ and dashed lines to fluctuations in ρ,

We represent in Fig.7 the spatio–temporal distribution of the energy density for scattering due to fluctuating μ (solid line) and for scattering due to fluctuating ρ (dashed line), which is obtained by integration of (25) and (27) over ϕ, For comparison, we have also depicted the spatio–temporal distribution of the energy density resulting from a single-scattering model based on the assumption of isotropic scatters (Sato 1977) (dotted line). In the three cases the scattering generates a logarithmic blow-up of the energy density near r=t, as is typical for single scattering (Sato & Fehler 1998). However, Fig.7(a) shows that the spatial distribution of coda energy is more uniform for the isotropic case than for the other cases. The energy density concentrated in the vicinity of the source for the scattering generated by the fluctuations of μ is about two times larger than for the others.

7

. (a) Spatial distribution at t=1 and (b) time evolution of energy density at r=1 for the single-scattering model. Solid lines correspond to media with fluctuations in Lamé parameter μ, dashed lines to media with fluctuations in mass density ρ and dotted lines to a media with isotropic scatterers.

The temporal evolution of the energy densities is depicted in Fig.7(b). In the inset, we show with a solid line (dashed line) the relative difference of the energy density between a medium with fluctuations only in μ (ρ) and the energy density of the single isotropic scattering model. It is clearly seen that the relative difference between the energy densities for these models is smaller in the early coda than in the late coda, where it grows monotonically. The relative difference asymptotically approaches 3/5 (1/3) for the first (second) case, since the asymptotic behaviour of I^μ(r,t) [I^ρ(r, t)] is 5/(4πt²) [3/(4πt²)], while the asymptotic behaviour for the single isotropic scatter-ing model is 1/(2πt²). Note that this is true independently of the source–receiver distance.

5 Numerical Results in Spherical Geometry

In this section we present Monte Carlo simulations for the energy density I and the degree of polarization α₁^{μ, ρ} when multiple scattering is taken into account. We show a detailed analysis of the full-space case. Some examples for the half-space and slab geometries are given in the next section.

In Fig.8 we plot the spatial energy density distribution at different times. As in the previous section, the solution for the medium with fluctuating Lamé parameter μ (first case) is represented with solid lines, and the solution for the medium with fluctuating density ρ (second case) with dashed lines. The dotted lines correspond to the scalar isotropic scattering model, which is shown for comparison. The difference between the three cases is apparent for short lapse times (Figs8a and b) and diminishes for larger times (Figs8c and d). The reason is that single scattering is dominant for short times and the spatial energy density distribution strongly depends on the underlying scattering law, as we showed in the previous section. When time increases, multiple scattering destroys this effect, and all curves tend to the same diffusion solution.

8

. Spatial energy density distributions at ×t=0.5, 1, 2 and 4 for a medium with fluctuations only in the Lamé parameter (solid lines) and for a medium with fluctuations only in the mass density (dashed lines). Comparison with isotropic scattering is shown by dotted lines.

More interesting for inversion problems (see e.g. Wu 1985; Fehler et al. 1992; Mayeda et al. 1992; Hoshiba 1995) are the time traces shown in Fig.9. Again, for small source–receiver distances, the early codas are not equal, since high-order multiple scattering is absent. However, the late codas look alike. For large hypocentral distances the whole time traces are similar. Hoshiba (1995) compared the early coda with the later one to estimate the amount of scattering in the forward half-space in western Japan. He found that coda wave energy concentrated just after the S-wave arrival and was more important than predicted by the isotropic scattering model, so asymmetry factors were introduced to explain this effect. Figs9(a)–(c) show that the coda envelope depends on the scattering phase matrix, and that more energy density is concentrated in the early coda for the isotropic scattering model than for the others. Consequently, to account for the large amount observed by Hoshiba (1995) in the early coda, we would also need to consider models that include forward scattering effects through power spectral densities R^_μμ and R^_ρρ that are no longer constant.

9

. Time traces at r=0.5, 1, 2 and 4 for a medium with fluctuations only in the Lamé parameter (solid lines) and for a medium with fluctuations only in the mass density (dashed lines). Comparison with isotropic scattering is shown by dotted lines.

Figs10 and 11 show the spatial distribution and time depen-dence of the degree of polarization, respectively. The solid lines represent case 1 and the dashed lines case 2. In Fig.10(a), we recognize the already known spatial distribution profiles analysed in Section 4 for single scattering. The degree of polarization for fluctuating μ shows three peaks, the second of them coinciding in position and magnitude with the unique one for fluctuating ρ, Nevertheless, it is seen in Figs10(b)–(d) that the closest peak to the source diminishes very quickly compared with the other two, and that the second peak is next to vanish. This effect is due to multiple scattering, which is responsible for depolarization. The closer we are to the source, the bigger the mean number of scattering events. On the other hand, the polarization due to the interaction between S waves and inhomogeneities in ρ (dashed lines) remains longer in the medium. This last fact can be easily observed in Fig.11, where we have plotted time traces for the degree of polarization. What is remarkable here, compared with Fig.6(b), is the fast decay of the polarization for the first case. This is explained in Fig.12, where we have represented the spatial distribution of the degree of polarization of particles having one, two, four and eight scattering events for both cases. This figure clearly shows that the depolarization for fluctuations in μ takes place much faster than the depolarization for fluctuations in ρ—only the waves closest to the first arrival remain more than 10 per cent polarized after eight scatterings. We point out that the discrepancy between the curve in Fig.12(a) for one scattering interaction and the analytical curve calculated in the single-scattering approximation (Fig.6a) is due to the discretization of the angular direction; that is, we have replaced the integral in eq.(28) by a sum over 12 discrete directions. We checked that the agreement between both figures became perfect by increasing the number of discrete directions.

10

. Same as Fig.8 but for the degree of polarization.

11

. Same as Fig.9 but for the degree of polarization.

12

. Spatial distribution of the degree of polarization of particles having one (solid line), two (dotted line), four (dash-dotted line) and eight (dashed line) scattering events for (a) R●_ρρ=0 and (b) R●_μμ=0.

Figs9 and 11 show that polarization is more sensitive to the scattering phase matrix arising from different forms of fluctuations than the energy density.

6 Free-surface effects on polarization

In this section we analyse the effect on the degree of polarization of a perfectly reflecting horizontal interface. The main question is by how much the degree of polarization is modified by the interaction between incident and reflected waves. We have addressed the main differences between the two scattering phase matrices coming from fluctuations in the Lamé parameter μ and the density ρ, We consider here only fluctuations in the density.

Fig.13 shows the degree of polarization for time t=2. The non-polarized isotropic point source is at a distance z=1 from the interface. Clearly, we can split the domain (r_{2− D}, z) into two areas. The lower part is actually identical to the spherical symmetry case, and corresponds to waves that did not see the interface. The maximum of the degree of polarization, about 0.5, perfectly agrees with the one shown in Fig.10(c). The upper part corresponds to the area where some energy has been reflected from the free surface. We observe that the degree of polarization is lower than in spherical geometry. This is because the polarization of the reflected part adds somewhat destructively to that of the part that did not see the interface.

13

. Degree of polarization in half-space geometry for a medium with fluctuations in the density. The non-polarized isotropic point source is at z=0 and the free surface is at z=1.

More interesting is the value of the degree of polarization at the surface, where it can be measured in practice. In Fig.14 it is seen that in spite of the relative diminution of polarization caused by the interface, there is still a significant amount of polarized energy, at least for early times. As in the spherical geometry, the degree of polarization starts from zero, then quickly grows until it reaches a maximum, and finally slowly decreases back to zero for large times.

14

. Degree of polarization at the surface for density fluctuations. The non-polarized isotropic point source is at a distance z=1 from the interface.

We have represented in Fig.15 the degree of polarization at the free surface in slab geometry (an infinite medium in the x–y directions, with finite thickness d in the z direction). In Fig. 15, d is equal to one mean free path, and the unpolarized source is 0.75 mean free paths away from the free surface (located at z=0). For simplification, the lower interface is also assumed to be perfectly reflecting. As in Fig. 14, a first peak appears approximately 0.25 mean free times after the first arrival reaches the zero-offset location. However, a second peak appears due to the energy that has been reflected from the lower boundary. Observe the presence of a minimum between these two maxima, where the degree of polarization is close to zero. This is explained by the unpolarized waves, which, travelling from the source to the lower boundary, are reflected and then received at the free surface.

15

. Degree of polarization pattern at the surface for the slab problem for a medium with fluctuations in the density. The non-polarized isotropic point source is at a distance 0.75 from the surface and 0.25 from the lower boundary.

7 Conclusions

We have studied the polarization and depolarization effects of S waves for media with fluctuating Lamé parameter μ or fluctuating mass density ρ, Coda energy envelopes are found to be very different for scattering caused by fluctuations in the density, scattering caused by fluctuations in the Lamé parameter and isotropic scattering, when single scattering is dominant compared to multiple scattering. Analytical solutions in the single-scattering approximation have been presented. We have underlined important polarization effects and demonstrated that the amount of energy density that is polarized after single scattering can be as large as 80 per cent, even for unpolarized sources. In addition, we have shown that these effects are quantitatively and qualitatively distinguishable for the two kinds of scattering phase matrices considered. Both the polarization and depolarization of S waves take place faster for fluctuating Lamé parameter than for fluctuating mass density. This suggests the use of polarization measurements of S waves to reconstruct the power spectrum of the fluctuations of the Lamé parameter μ and the mass density ρ in the lithosphere from real seismic data. We refer to McCormick (1990) for the reconstruction of properties of the medium from boundary measurements including polarization. We have also presented a Monte Carlo method for S waves that is easy to implement and capable of addressing general 3-D geometries. For half-space and slab problems, we have shown that the presence of an interface lowers the degree of polarization. However, this degree remains clearly observable.

Acknowledgements

We would like to thank Joseph B. Keller, George Papanicolaou and Leonid Ryzhik for interesting discussions. We acknowledge helpful comments and criticisms from the anonymous referees, which helped us to improve the readability of the original manuscript. This work was supported by AFOSR grant 49620-98-1-0211 and NSF grant 9709320.

Appendices

Appendix A: Probability Scattering density

We give here an explicit expression for eqs (18) and (19) in the case of fluctuations in the density. Let S=T in these expressions. The components of the T matrix (8) are given by

Upon expanding the trace of TWT^*, we have

where C=3/(8π). From eq.(19) we have

For S=Γ, corresponding to fluctuations in the Lamé parameter μ, expressions (A2) and (A3) are obtained by substituting T_ij for Γ_ij [(i, j)=1,2] given by (8), and C=5/(8π).

Appendix B: Energy Density and Degree of Polarization in the Single-Scattering Approximation

The energy density that arrives at a point r after scattering at r₁ due to inhomogeneities in ρ (see Fig.2) is given by

The factors 1/4πr₁² and 1/4πr₂² in (B1) come from the geometrical spreading of the energy density from the source to the scatter and from this to the receiver, respectively. To obtain the total energy density at time t we have to integrate over the whole space under the condition that the total traveltime of the wave is equal to t. Then,

We note that

Therefore, we have

where H is the step function. From eq.(B3) we infer the energy density arriving with angle ϕ given in eq.(27). The numerator in (29) is the same as (B3) except that the positive sign in front of cos²Θ is replaced by a negative sign.

To integrate I(r, t) analytically it is more convenient to change to prolate spheroidal coordinates (νω, φ) directly in (B2) as in Sato & Fehler (1998). The change of variables is as follows:

where r is the distance between the source and the receiver, and ω∈[− 1, 1], ν∈[1, ∞] and φ∈[0, 2π]. The Jacobian of this transformation is (rr₁r₂)/2dνdωdφ and

With these relations we end up with the integration of a rational function, which can be performed trivially. The result for the total energy density and the degree of polarization in the case of fluctuating ρ is

where ν=t/r (ν=vt/r in dimensional units). The total energy density for inhomogeneities in μ can also be calculated by this procedure. The result is

However, the integral for the degree of polarization α₁^μ(r, t) is more complicated because (24) changes sign over the interval of integration. For inhomogeneities in μ, we have evaluated eq.(29) numerically.

Appendix C: Implementation of the Monte Carlo Method

In this section, we describe the implementation of the Monte Carlo method for solving the radiative transfer equation (13).

C1 Decomposition in the number of scattering events

Let us denote by

the probability that a particle starting at time 0 from z₀=(x₀, k₀, p₁⁰, p₂⁰) be at time T at z=(x, k, p₁, p₂)about dz=dxdkdp₁dp₂. We can decompose P(T, z, dz;z₀) over the number of scattering events between 0 and T as

where N_s is the number of scattering events, or jumps, of the process (X(t),K(t), P(t)) between 0 and T. Here | means conditional expectation and P(T, N_s=n) is the probability of having exactly n scattering events between 0 and T. In other words, the total probability density is decomposed over the number of scattering events that the process experiences between 0 and T. The Poisson distribution of the law of scattering times (17) allows one to compute each term in the sum (C2) separately by the Monte Carlo method. We can check that

The times at which the n scatterings occur are found by drawing n values uniformly between 0 and T and arranging them in increasing order.

The first advantage in computing each term of the sum successively is the reduction of the statistical fluctuations. Indeed, we include some information on the exact proportion of single scattering, double scattering, etc., in the numerical solution. Moreover, we have access to these quantities and can easily find out in which areas single scattering or multiple scattering is dominant. The second advantage is architectural. Our code is written in _Matlab, which is much faster when properly vectorized. We calculate in turn the contributions of particles having exactly n scattering events during a fixed time T. Here we have 1 ≤n≤K, where K is chosen so that the particles having more than K interactions during an interval of time of size T are negligible, i.e.

is small. Once n is fixed, the number of operations required to compute the trajectories of M particles is independent of the particle and therefore easily vectorized. The total number of particles N=N_loops×M is computed by N_loops successive calculations of packets of M particles.

For a fixed n, the operations required to simulate each trajectory are as follows. We first compute the times at which the process scatters. For each scattering we determine a new velocity k and new polarization parameters q and u of the outgoing particle according to (18) and (19). As seen in (18), the probability distribution of the outgoing velocity depends on the incoming velocity, on the incoming polarization parameters and on the outgoing velocity. It is, however, independent of the outgoing polarization parameters. The outgoing polarization parameters are determined by (19) once the outgoing velocity is known. Between scattering events, particles travel along straight characteristics with velocity v=1.

C2 Spatial and angular discretization

Due to the large number of computations involved tosimulate each trajectory, we have discretized the velocity and polarization spaces (144 discrete velocities and 16 discrete polarization states). We precompute the probability distribution of the outgoing parameters for these discrete values of the velocity and polarization. We then replace the exact probability distribution for the outgoing velocity by a piece-wise constant probability distribution. Each interaction of a particle with the underlying medium is sampled accordingly.

At final time, the particles must be boxed to perform the averaging over the polarization parameters. We discretize the spatial domain with n_x=50^d boxes B_i, 1 ≤i≤n_x, and the angular domain with n_k=144 directions k_l, 1 ≤l≤n_k. Here d is the effective dimension. In spherical geometry, the effective dimension d=1 due to spherical invariance. For the slab problem, the effective geometry is d=2 due to cylindrical invariance. Note that to use the geometrical symmetries, the polarization parameters have to be rotated properly (we refer to Chandrasekhar 1960 for the details). Each particle n for 1 ≤n≤N is represented at final time T by (x(n), k(n), q(n), u(n)), which are its position, velocity and polarization parameters, respectively. Then, for each i and l we compute the quantities

where

The quantities of interest are the energy density and the degree of polarization defined for each box by

The number of discrete velocities, discrete polarization states and spatial boxes has been tailored so as to have perfect visual agreement between our simulations and the analytical results obtained for the single-scattering approximation and in the long time diffusion regime. Finer discretizations in specific configurations have revealed no substantial improvement. The number of generated trajectories has been chosen so as to obtain small statistical fluctuations in every box of the discretization.

We thus believe that our numerical simulations are close to the exact solutions of the radiative transfer equation (13).