TThe combinatorics of scattering in layered media
Peter C. Gibson ∗ May 17, 2013
Abstract
Reflection and transmission of waves in piecewise constant layeredmedia are important in various imaging modalities and have been stud-ied extensively. Despite this, no exact time domain formulas for theGreen’s functions have been established. Indeed, there is an under-lying combinatorial obstacle: the analysis of scattering sequences. Inthe present paper we exploit a representation of scattering sequencesin terms of trees to solve completely the inherent combinatorial prob-lem, and thereby derive new, explicit formulas for the reflection andtransmission Green’s functions.
The analysis of wave propagation in layered media is a well-developed sub-ject, with applications to acoustic, electromagnetic and seismic imaging; see[2], [4], [1] and the numerous references therein. In the present introductionwe describe briefly the physical framework and cite some needed facts. Inthe next section we explain a fundamental combinatorial problem arisingfrom the given framework—the enumeration of scattering sequences. Themain contribution of the present paper is to solve the combinatorial prob-lem completely, thereby producing new, exact formulas for the reflectionand transmission Green’s functions. Our main results are Theorem 2.1 andTheorem 3.1.The basic framework is as follows. Let ( x, y, z ) denote euclidean coor-dinates for a solid three-dimensional acoustic medium whose physical pa-rameters (density and bulk modulus) vary only in the z -direction. Suppose ∗ Dept. of Mathematics & Statistics, York University, 4700 Keele St., Toronto, Ontario,Canada, M3J 1P3, pcgibson @ yorku . ca a r X i v : . [ m a t h . C O ] M a y eter C. Gibson The combinatorics of scattering z , hav-ing jumps at the M + 1 locations z < z < · · · < z M . Thus the solid contains M homogenous layers, the n th layer correspondingto the z -interval z n − < z < z n ;the layers are sandwiched between the two half spaces z < z and z >z M . We shall refer to the z -direction as depth, and depict it as increasingdownward, as in Figure 1. z (cid:45) z (cid:45) z z z z z z z z z z z z (cid:187)(cid:187) z M (cid:45) z M (cid:45) z M z M z M (cid:43) z M (cid:43) Source PulseSource Pulse Reflected SignalReflected SignalTransmitted SignalTransmitted Signal
Figure 1: A layered medium with M layers, sit-ting between two half spaces. The depth coordi-nate z increases downward.From the perspective of acoustic imaging the problem is to infer thephysical parameters for the M layers by probing them with acoustic wavesin one of two ways. The first way is to send an acoustic pulse from a fixeddepth z − < z toward the interface at z and to record the pulse trainreflected back from the M layers as it crosses the original depth z − . The reflection problem is to infer the structure of the M layers from this recordedreflection data. A second experiment is to transmit a pulse from z − asbefore, but to record the resulting pulse train that is transmitted across the M layers to some fixed depth z M +1 > z M . The transmission problem is to17 May 2013 eter C. Gibson The combinatorics of scattering M layers from this transmitted data, recorded at z M +1 . It is assumed for both the reflection and transmission experimentsthat the initial pulse is a plane wave of the form f ( z − ct ), where t denotestime and c ( z ) denotes the speed of sound at depth z . Thus although thephysical setting is three-dimensional, by symmetry the analysis reduces toone spatial dimension, the z -direction.If the initial pulse is idealized to a Dirac delta function initially centredat z − , so of the form δ ( z − z − − ct ) , then the reflection data G ( t ) and the transmission data H ( t ) are Green’sfunctions for the governing partial differential equation (see Appendix A),and they have the form G ( t ) = ∞ (cid:88) n =1 a n δ ( t − σ n ) , (1.1) H ( t ) = ∞ (cid:88) n =1 b n δ ( t − σ n ) . (1.2)The numbers a n and b n will be referred to as reflection amplitudes and trans-mission amplitudes, respectively; the numbers σ n and σ (cid:48) n will be referred toas arrival times.As a consequence of the standard theory, both G and H are completelydetermined by two vectors:1. the sequence R = ( R , . . . , R M ) of reflection coefficients at the layerboundaries; and2. the sequence τ (cid:48) = ( τ , . . . , τ M +1 ) of two-way travel times, where τ n denotes twice the time it takes a downward-traveling acoustic wave togo from z n − to z n .The travel time τ M +1 is obviously irrelevant for the reflection problem, so G is in fact determined by R and τ = ( τ , . . . , τ M ). We incorporate this basicfact into our notation, writing G ( τ,R ) and H ( τ (cid:48) ,R ) to denote the respective reflection and transmission data determined by agiven pair ( τ (cid:48) , R ).Different layered media that correspond to the same parameters ( τ (cid:48) , R )are for present purposes indistinguishable, and we shall simply refer to the17 May 2013 eter C. Gibson The combinatorics of scattering τ (cid:48) , R ) as a medium, letting it be understood that a class of physicalsystems are thereby represented. (See Appendix A.)The above facts lead naturally to a basic question: What is the formulafor the amplitude coefficients a n and b n in terms of ( τ (cid:48) , R )? It turns outthat no finite, closed form formula has ever been established—only seriesexpansions that must be estimated when it comes to practical computation.(See, for example, [4, Chapter 3], [1, Chapter 2], [2, Section 2.5]; concerningtheoretical developments applied specifically to seismic, see [8], [9], [11],[5])This is because there is a substantial combinatorial problem blocking theway to an exact formula. The purpose of the present paper is to solve thecombinatorial problem directly, and to present the resulting explicit formulasfor the reflection and transmission Green’s functions G and H in terms of( τ (cid:48) , R ). The one-dimensional wave equation for a homogenous medium generatestraveling wave solutions; in a layered medium one has to take into accountthe behaviour of these traveling waves at interfaces between homogenouslayers. This behaviour depends on the reflection coefficient R n associatedwith the interface at z n in the following way. When a wave f ( z − z n − c n t )traveling downward from z n − toward z n for t < t = 0,it splits into a reflected wave, R n f ( z − z n + c n t )that travels back up toward z n − , and a transmitted wave (cid:112) − R n f ( z − z n − c n +1 t )that continues down toward z n +1 (at modified speed c n +1 ). The reflectedand transmitted waves then hit the respective interfaces at z n − and z n +1 ,generating two new sets of reflected and transmitted waves, and a cascade ofsuccessive reflections and transmissions continues indefinitely. The case of aninitial wave traveling upward from z n +1 toward z n instead of downward from z n − is similar, except that the reflection coefficient − R n applies instead of R n ; the transmission coefficient (cid:112) − R n remains unchanged. We use thenotation T = ( T , . . . , T M ) for the transmission coefficients that correspondto reflection coefficients R = ( R , . . . , R M ) by way of the formula T n = (cid:112) − R n (0 ≤ n ≤ M ) .
17 May 2013 eter C. Gibson The combinatorics of scattering
5A standard idea is to view G ( τ,R ) ( t ) as the sum total of all possiblesequences of successive reflections and transmissions of an initial pulse at z − that eventually return to z − . (See [4, Chapter 3] for details.) Theidea is worth illustrating further, since it underpins the arguments below.For example, consider an initial downward traveling unit pulse of the form δ ( z − z − − c t ). After τ / z . Part of it is transmitted into the first layer as T δ (cid:0) z − z − c ( t − τ ) (cid:1) ,which, after another τ / z . Part of this pulseis reflected back into first layer as R T δ (cid:0) z − z + c ( t − τ + τ ) (cid:1) . Travelingback up to z , the latter reaches z after a further τ / z − as R T δ (cid:0) z − z + c ( t − τ +2 τ ) (cid:1) , (2.1)arriving at z − at time σ = τ + τ . Thus the part of the initial pulse thattraverses the sequence of depths p = ( z − , z , z , z , z − ) returns to z − withan amplitude α = R T at arrival time σ = τ + τ , thereby contributing aterm of the form αδ ( t − σ ) to G ( τ,R ) ( t ). The sequence p is called a scatteringsequence, and the amplitude α is called the weight of p . As mentionedearlier, the impulse response itself is a delta train of the form G ( τ,R ) ( t ) = ∞ (cid:88) n =1 α n δ ( t − σ n ) , (2.2)composed of the cumulative contributions of all possible scattering sequencesreturning to z − , with their associated weights and arrival times. In theabove example, p is the only scattering sequence having arrival time σ . Butin general different scattering sequences may arrive simultaneously, so eachamplitude α n occurring in (2.2) is the sum of the weights of all scatter-ing sequences arriving at time σ n . Given an arrival time σ n , the essentialcombinatorial problem, then, is to enumerate all the associated scatteringsequences, along with their weights, to derive a formula for α n . The first step is to establish some terminology and notation, as follows. Ascattering sequence that starts and ends at z − is represented by a path inthe graph z (cid:45) z (cid:45) z z z z z z z z z z ...... z M (cid:45) z M (cid:45) z M z M (2.3)17 May 2013 eter C. Gibson The combinatorics of scattering M ≥
1, let S M denote the set of all sequences ofthe form p = ( p , p , . . . , p L )such that L ≥ p = p L = z − and p n ∈ { z , . . . , z M } if 1 ≤ n ≤ L −
1; (2.4a) ∀ n with 0 ≤ n ≤ L − , ∃ j with − ≤ j ≤ M − { p n , p n +1 } = { z j , z j +1 } . (2.4b)The elements of S M will be referred to as scattering sequences. Condition(2.4a) says a scattering sequence starts and ends at z − , and condition (2.4b)(which refers to un ordered pairs) says that adjacent terms in a scatteringsequence are adjacent vertices in the graph (2.3).For example, for 0 ≤ n ≤ M , the shortest scattering sequence thatreaches z n is ( z − , z , z , . . . , z n − , z n , z n − , . . . , z , z , z − ); (2.5)this is called a primary scattering sequence. A scattering sequence reachingmaximum depth z n that is not shortest possible is called a multiple scatteringsequence. The weight w ( p ) corresponding to a scattering sequence p = ( p , p , . . . , p L )in an M -layer medium ( τ, R ) is defined as follows. For each n in the range1 ≤ n ≤ L −
1, and given that p n = z j , define w n = R j if p n − = p n +1 = z j − − R j if p n − = p n +1 = z j +1 (cid:113) − R j otherwise . (2.6)The three possibilities correspond respectively to: reflection inside the j thlayer at z j ; reflection inside the ( j +1)st layer at z j ; and transmission betweenthe j th and ( j + 1)st layers. Finally, set w ( p ) = L − (cid:89) n =1 w n . Thus the part of an initial unit impulse that traverses p returns to z − withamplitude w ( p ).17 May 2013 eter C. Gibson The combinatorics of scattering We define two maps, κ, β : S M → Z M +1 that associate integer vectors to a given scattering sequence.A scattering sequence p ∈ S M may be represented graphically as inFigure 2; Stanley [10] calls such a representation a Dyck path. z (cid:45) z (cid:45) z z z z z z z z z z z z z z z z z z Figure 2: The Dyck path for a scattering sequence p , with the horizontalcoordinate now representing time, which increases to the right. z (cid:45) z (cid:45) z z z z z z z z z z z z z z z z z z Figure 3: The intervals I j , in red. There are five intervals, two of whichhave positive length, so k = 5 and b = 2.Given the Dyck path for a scattering sequence p ∈ S M , let t denote thehorizontal coordinate and (as usual) let z denote the vertical coordinate.Let U denote the portion of the t, z -plane on or above the Dyck path and at17 May 2013 eter C. Gibson The combinatorics of scattering z = z − —the shaded region in Figure 3. For each n in the range0 ≤ n ≤ M , consider the horizontal line L n at depth z n . The intersection L n ∩ U consists of a disjoint union of closed intervals I nj , where 1 ≤ j ≤ k n ;see Figure 3. The intervals of I nj are of two types: degenerate intervalsconsisting of a single points; and non-degenerate intervals having positivelength. Letting k n denote the total number of intervals I nj , and letting b n denote the number of non-degenerate intervals, set κ ( p ) = k = ( k , . . . , k M ) and β ( p ) = b = ( b , . . . , b M ) . Observe that the entry k n of the vector k = κ ( p ) counts the number oftimes the Dyck path crosses back and forth across the n th layer z n − < z A tree is a connected cycle-free graph. The vertices of a tree are divided intothree anatomical types, as follows: the root is a single, specially designatedvertex; a non-root that belongs to just one edge is called a leaf; all othernon-root vertices are called branch points. Vertices in a tree have a height,determined by their distance (in the sense of shortest path) to the root. SeeFigure 4. z (cid:45) z (cid:45) z z z z z z z z z z z z z z z z z z Figure 4: A tree, with the leaves coloured green. The root is atthe bottom, and horizontal lines indicate the various heightsof vertices.The association between a scattering sequence and a tree is well-known(see [10, Exercise 6.19]) and arises, for instance, in the analysis of Brownianexcursions and superprocesses (see [7, Section 1.1]). The tree representing ascattering sequence p ∈ S M may be obtained simply by collapsing its Dyckpath, as follows. Recall the intervals I nj used to define κ ( p ) and β ( p ) above;for present purposes let I − denote the intersection of z = z − with U . Tocollapse the Dyck path, contract each of the intervals I nj to a point, keepingthe distances between intervals unchanged, and interpolate this horizontalcontraction linearly on each depth interval z n − < z < z n . This opera-tion transforms the original Dyck path into a tree (in fact it is an isotopybetween the region U and the resulting tree). See Figure 5. Note that the17 May 2013 eter C. Gibson The combinatorics of scattering I nj (coloured green for emphasis) end up as leaves, whilenon-degenerate intervals are contracted to branch points of the tree—exceptfor I − , which is contracted to the root.Conversely, given a tree, one may recover the original scattering sequenceby tracing the outline of the tree, from the root (keeping the tree on theleft), and recording the depths of the vertices in the order that they arepassed.(a) z (cid:45) z (cid:45) z z z z z z z z z z z z z z z z z z (b) z (cid:45) z (cid:45) z z z z z z z z z z z z z z z z z z (c) z (cid:45) z (cid:45) z z z z z z z z z z z z z z z z z z (d) z (cid:45) z (cid:45) z z z z z z z z z z z z z z z z z z (e) z (cid:45) z (cid:45) z z z z z z z z z z z z z z z z z z (f) z (cid:45) z (cid:45) z z z z z z z z z z z z z z z z z z Figure 5: The collapsing of a Dyck path (a) to a tree (f) (depicted upsidedown), an operation which is reversible.Observe that the vectors ( k, b ) = ( κ ( p ) , β ( p )) have a simple interpreta-tion in terms of the tree representing p ∈ S M . For each 0 ≤ n ≤ M , k n isthe number of vertices at depth z n , and b n is the number of branch pointsat z n . (This is the reason for calling b the branch count vector for p .) Thereare some evident constraints on these quantities. Note first that b M = 0.Furthermore, for 0 ≤ n ≤ M − { , k n +1 } ≤ b n ≤ k n +1 (0 ≤ n ≤ M − , (2.10)since each vertex at z n +1 is connected by an edge to a unique branch pointat z n . It is convenient to refer to the left shift ˜ k of k = κ ( p ); that is,˜ k = ( k , k , . . . , k M , ∈ Z M +1 . (2.11)17 May 2013 eter C. Gibson The combinatorics of scattering { , ˜ k n } ≤ b n ≤ min (cid:8) k n , ˜ k n (cid:9) (0 ≤ n ≤ M ) . (2.12)An easy application of the tree representation is to determine the possiblevalues of b = β ( p ) for each given k = κ ( p ), or in other words, to determinethe set V ( k ) = β ( κ − ( k )) . (2.13)In fact V ( k ) is determined precisely by (2.12). Proposition 2.1 Given k ∈ L M , V ( k ) = (cid:8) b (cid:12)(cid:12) min { , ˜ k } ≤ b ≤ min { k, ˜ k } (cid:9) ; equivalently, V ( k ) may be expressed as a Cartesian product of sets, V ( k ) = V × V × · · · × V M , where V n = (cid:8) min { , ˜ k n } , min { , ˜ k n } + 1 , . . . , min { k n , ˜ k n } (cid:9) (0 ≤ n ≤ M ) . Here ∈ Z M +1 is the vector whose entries are all 1. The minimum is to beinterpreted entrywise, meaning that for x, y ∈ R M +1 ,min { x, y } = (cid:0) min { x , y } , min { x , y } , . . . , min { x M , y M } (cid:1) ∈ R M +1 . Proof . The constraints (2.12) imply that V ( k ) ⊂ V ×· · ·× V M . It remains toshow that any b ∈ V × · · · × V M belongs to V ( k ), which entails showing thatthere exists a scattering sequence p such that ( κ ( p ) , β ( p )) = ( k, b ). It sufficesto construct a tree representing such a p , as follows. Given b ∈ V ×· · ·× V M ,place k n vertices, consisting of b n branch points and k n − b n leaves (in anyorder), at depth z n , for 0 ≤ n ≤ M , and place a root at z − (consideredas a branch point). Let N denote the largest index such that k N (cid:54) = 0. For0 ≤ n ≤ N , 1 ≤ b n − ≤ min { k n − , k n } , so there exists a surjection f n from the k n vertices at z n onto the b n − branchpoints at z n − . The tree representing p is completed by drawing an edgefrom each vertex v at z n to f n ( v ).17 May 2013 eter C. Gibson The combinatorics of scattering The tree representation facilitates deriving a simple formula for the weights.The following lemma uses multi-index notation, whereby given a vector s =( s , s , . . . , s M ) and an integer vector d = ( d , . . . , d M ), s d = M (cid:89) n =0 s d n n . Lemma 2.2 Let p ∈ S M be a scattering sequence in an M -layer medium ( τ, R ) , and set ( k, b ) = ( κ ( p ) , β ( p )) . Then w ( p ) = ( − R ) ˜ k − b R k − b T b . Proof . Consider the tree representing p . Observe that each instance of w n = R j in (2.6) corresponds to a unique leaf at z j (see Figure 5). Sincethere are k n − b n leaves at z n this results in a total contribution of R k − b .Let v be a branch point at z j having d v edges to vertices at z j +1 . Observethat precisely d v − w n = − R j in (2.6), and every occurrence of w n = − R j arises this way. The sum totalof numbers d v − v at z j is simply k n +1 − b n = ˜ k n − b n ,making for a total contribution over all depths of ( − R ) ˜ k − b .Finally, each instance of transmission from the j th layer to the ( j + 1)stlayer in p corresponds to a vertex in the tree representing p which is nota leaf, i.e., to a branch point at z j ; and, since the path p starts and endsat z − every such transmission has a corresponding return transmission inthe opposite direction, from the ( j + 1)st layer to the j th layer. There are b j branch points at z j , and each of these corresponds to two transmissionsacross the boundary at z j , making for a total contribution to w ( p ) of T b .Since every w n is covered by one of the above cases, the lemma follows. A formula for the coefficients a ( R, k ) in (2.8), defined as the summation(2.9), is now within easy reach. Recall that the binomial coefficient (cid:0) xy (cid:1) fora pair of non-negative integer vectors x, y ∈ Z M +1 , with y ≤ x , is to beinterpreted as (cid:18) xy (cid:19) = M (cid:89) n =0 (cid:18) x n y n (cid:19) . (The inequality y ≤ x means that x − y has non-negative entries.)17 May 2013 eter C. Gibson The combinatorics of scattering Lemma 2.3 Let k ∈ L M and let min { , ˜ k } ≤ b ≤ min { k, ˜ k } . Then { p ∈ S M | ( κ ( p ) , β ( p )) = ( k, b ) } = (cid:18) kb (cid:19)(cid:18) ˜ k − ub − u (cid:19) . Proof . To count the number of scattering sequences having a given transitcount vector, it suffices to count the number of corresponding trees—whichis straightforward. Consider first the arrangement of vertices in a tree forwhich ( κ ( p ) , β ( p )) = ( k, b ). At each depth z n , there are k n vertices of which b n are branch points and k n − b n are leaves. There are (cid:0) k n b n (cid:1) ways of arrangingthese from left to right, making for a total of (cid:18) kb (cid:19) = M (cid:89) j =0 (cid:18) k n b n (cid:19) (2.14)possible vertex arrangements. (There is only one way to place the root at z − , which may be ignored.)For each vertex arrangement there are various possible edge arrange-ments, as follows. Each of the k n +1 = ˜ k n vertices at z n +1 must be con-nected by an edge to one of the b n branch points at z n , respecting thevertex ordering (so that edges don’t cross). This is equivalent to choosinga b n -part ordered partition of the integer ˜ k n . If ˜ k n ≥ 1, there are (cid:0) ˜ k n − b n − (cid:1) possible choices; and if ˜ k n = 0 then b n = 0 and there is 1 = (cid:0) ˜ k n b n (cid:1) (empty)arrangement. Letting N denote the largest index for which ˜ k N (cid:54) = 0, thetotal number of edge arrangements is (cid:18) ˜ k − ub − u (cid:19) = N (cid:89) n =0 (cid:18) ˜ k n − b n − (cid:19) . (2.15)Combining (2.14) and (2.15) yields a total tree count of { p ∈ S M | ( κ ( p ) , β ( p )) = ( k, b ) } = (cid:18) kb (cid:19)(cid:18) ˜ k − ub − u (cid:19) , completing the proof.Combining the foregoing lemmas gives a formula for the reflection Green’sfunction, as follows. Theorem 2.1 Let ( τ, R ) be an M -layer medium for some integer M ≥ .Then G ( τ,R ) ( t ) = (cid:88) k ∈ L M a ( R, k ) δ ( t − (cid:104) k, τ (cid:105) )17 May 2013 eter C. Gibson The combinatorics of scattering where for each k ∈ L M , the amplitude a ( R, k ) is given by the followingformula. Setting u = min { , ˜ k } , a ( R, k ) = (cid:88) b ∈ V ( k ) (cid:18) kb (cid:19)(cid:18) ˜ k − ub − u (cid:19) ( − R ) ˜ k − b R k − b T b , (2.16) where V ( k ) denotes the set of b ∈ Z M +1 such that u ≤ b ≤ min { k, ˜ k } .Proof . The total amplitude resulting from scattering sequences having agiven transit count vector k ∈ L M is a ( R, k ) = (cid:88) p ∈ S M such that κ ( p )= k w ( p ) . By Proposition 2.1 and Lemma 2.2 the above sum may be rearranged as a ( R, k ) = (cid:88) b ∈ V ( k ) (cid:88) p ∈ S M such that( κ ( p ) ,β ( p ))=( k,b ) w ( p )= (cid:88) b ∈ V ( k ) (cid:16) { p ∈ S M | ( κ ( p ) , β ( p )) = ( k, b ) } (cid:17) ( − R ) ˜ k − b R k − b T b . Applying Lemma 2.3 then gives the stated formula.Note that since each transmission coefficient T n = (cid:112) − R n occurs to aneven power in (2.16), the amplitude a ( R, k ) is a polynomial in the variables R , . . . , R M of precise degree (cid:104) k + ˜ k, (cid:105) . Note further that this polynomial depends only on the R n correspondingto the support of the vector k : if k n = 0, then a ( R, k ) does not dependon R n . Because the amplitudes are polynomial functions of the reflectioncoefficients, Theorem 2.1 is straightforward to code. This makes it possibleto compute G ( τ,R ) ( t ) exactly up to some cutoff time T > 0, and to do sovery efficiently—see Section 4.The simplest possible transit count vectors are those that correspond toprimary scattering sequences (2.5); given n ≥ 1, write k n for the primarytransit count vector defined as k nj = (cid:26) ≤ j ≤ n n < j . 17 May 2013 eter C. Gibson The combinatorics of scattering V ( k n ) = { k n − } is a singleton,so that the general formula (2.16) reduces just to a ( R, k n ) = ( − R ) R k n − k n − T k n − = R n n − (cid:89) j =0 (1 − R j ) , which is exactly Kunetz’s formula. The case of transmission is in many ways similar to that of reflection. Wetherefore give a much more compressed presentation, emphasizing only thosepoints where there is a substantial difference. In the case of transmission, a scattering sequence is a path in the graph z (cid:45) z (cid:45) z z z z z z z z z z ...... z M (cid:45) z M (cid:45) z M z M z M (cid:43) z M (cid:43) (3.1)that starts at z − and ends at z M +1 . (See Figure 1.) Formally, given aninteger M ≥ 1, let S (cid:48) M denote the set of all sequences of the form p = ( p , p , . . . , p L )such that: p = z − , p L = z M +1 and p n ∈ { z , . . . , z M } if 1 ≤ n ≤ L − 1; (3.2a) ∀ n with 0 ≤ n ≤ L − , ∃ j with − ≤ j ≤ M such that { p n , p n +1 } = { z j , z j +1 } . (3.2b)The elements of S (cid:48) M will be referred to as scattering sequences (or transmis-sion scattering sequences for emphasis).The weight of a scattering sequence in S (cid:48) M is computed exactly as for ascattering sequence in S M ; see Section 2.1.1.17 May 2013 eter C. Gibson The combinatorics of scattering z (cid:45) z (cid:45) z z z z (cid:187)(cid:187) z M (cid:45) z M (cid:45) z M z M z M (cid:43) z M (cid:43) Figure 6: The scattering path for a transmission scattering sequence. Aswith a Dyck path, the horizontal coordinate represents time, increasing tothe right. Figure 6 depicts what we call a scattering path, the analogue of a Dyck pathfor a transmission scattering sequence.Note that there are subpaths of a transmission scattering path that them-selves are Dyck paths. We call these Dyck subpaths. Collapsing the Dycksubpaths turns the scattering path into a tree, which has some additionalstructure. There is a distinguished leaf, the single leaf of maximum height(corresponding to z M +1 ), which we call the tip . The subpath leading directlyfrom the root to the tip will be referred to as the trunk . See Figure 7.Collapsing of Dyck paths is reversible, so there is a one-to-one correspon-dence between scattering paths and trees; the latter serve as a convenientrepresentation for the purpose of counting. Since every tree that correspondsto a scattering sequence has a trunk (including the root and the tip of thetree), it simplifies the combinatorial analysis to focus only on the remainingpart of the tree. As in Section 2, we define two maps, κ (cid:48) , β (cid:48) : S (cid:48) M → Z M +1 that associate integer vectors to a given transmission scattering sequence17 May 2013 eter C. Gibson The combinatorics of scattering z (cid:45) z (cid:45) z z z z (cid:187)(cid:187) z M (cid:45) z M (cid:45) z M z M z M (cid:43) z M (cid:43) z (cid:45) z (cid:45) z z z z (cid:187)(cid:187) z M (cid:45) z M (cid:45) z M z M z M (cid:43) z M (cid:43) Figure 7: The scattering path with Dyck subpaths filled in, and it’s as-sociated tree (depicted upside down, with the root at the top) obtained bycollapsing the Dyck subpaths. The tip is coloured yellow; the other leavesare green. p ∈ S (cid:48) M . For each n in the range 0 ≤ n ≤ M , define k n to be the numberof nodes at height z n not on the trunk. (Since there is precisely one nodeat each level belonging to the trunk, k n is one less than the total number ofnodes at height z n .) Set κ (cid:48) ( p ) = ( k , k , . . . , k M ) . Note that with this defintion k = 0; the other k n can take any non-negativeinteger value. Define m n to be the number of branch points at z n not onthe trunk, and set β (cid:48) ( p ) = ( m , . . . , m M ) . As before, we call k = ( k , . . . , k M ) a transit count vector and m = ( m , . . . , m M )a branch count vector, and we define ˜ k to be the left shift of k :˜ k = ( k , k , . . . , k M , . Observe that for each n , 0 ≤ m n ≤ min { k n , ˜ k n } . (3.3)Furthermore, letting Z + denote the nonnegative integers, the following resultis straightforward (we write 0 for both the number and the zero vector). Proposition 3.1 Given any k ∈ { } × Z M + and any ≤ m ≤ min { k, ˜ k } ,there exists a transmission scattering sequence p ∈ S (cid:48) M such that (cid:0) κ (cid:48) ( p ) , β (cid:48) ( p ) (cid:1) = ( k, m ) . 17 May 2013 eter C. Gibson The combinatorics of scattering Lemma 3.2 Let p ∈ S (cid:48) M be a transmission scattering sequence in an M -layer medium ( τ (cid:48) , R ) , and set ( k, m ) = ( κ (cid:48) ( p ) , β (cid:48) ( p )) . Then w ( p ) = ( − R ) ˜ k − m R k − m T m + . Proof . Consider the tree representing p . Observe that each instance of w n = R j in (2.6) corresponds to a unique leaf at z j (see the right-hand partof Figure 7). Since there are k n − m n leaves at z n this results in a totalcontribution of R k − m .Let v be a branch point at z j (possibly on the trunk) having d v edgesto vertices at z j +1 . Observe that precisely d v − w n = − R j in (2.6), and every occurrence of w n = − R j arises this way. The sum total of numbers d v − v at z j is simply ( k n +1 + 1) − ( m n + 1) = ˜ k n − m n , making for a total contributionover all depths of ( − R ) ˜ k − m .Finally, each instance of transmission from the j th layer to the ( j + 1)stlayer in p corresponds to a vertex in the tree representing p which is nota leaf, i.e., to a branch point at z j . Each branch point not on the trunkcorresponds to two transmissions across the boundary at z j (as for a Dyckpath), making for a total contribution to w ( p ) of T m . Each branch point onthe trunk corresponds to a single upward transmission, making for a totalcontribution of T . Since every w n is covered by one of the above cases, thelemma follows. The analogue of Lemma 2.3 for transmission scattering sequences is slightlysimpler, as follows. Lemma 3.3 Let k ∈ { } × Z M + and let ≤ m ≤ min { k, ˜ k } . Then (cid:8) p ∈ S (cid:48) M | ( κ (cid:48) ( p ) , β (cid:48) ( p )) = ( k, m ) (cid:9) = (cid:18) km (cid:19)(cid:18) ˜ km (cid:19) . Proof . To count the number of scattering sequences having a given tran-sit count vector, it suffices to count the number of corresponding trees,as before. Consider first the arrangement of vertices in a tree for which( κ (cid:48) ( p ) , β (cid:48) ( p )) = ( k, m ). To begin with, there is the trunk, which has fixedstructure. At each depth z n , there are k n vertices off the trunk, of which m n 17 May 2013 eter C. Gibson The combinatorics of scattering k n − m n are leaves. There are (cid:0) k n m n (cid:1) ways of arrangingthese from left to right, making for a total of (cid:18) km (cid:19) = M (cid:89) j =0 (cid:18) k n m n (cid:19) (3.4)possible vertex arrangements.For each vertex arrangement there are various possible edge arrange-ments, as follows. Each of the k n +1 = ˜ k n vertices at z n +1 must be connectedby an edge to one of the m n + 1 branch points at z n (including the branchpoint on the trunk), respecting the vertex ordering (so that edges don’tcross). This is equivalent to choosing a m n + 1-part ordered partition of theinteger ˜ k n , where the last part—corresponding to the trunk branch point—can be empty. There are (cid:0) ˜ k n m n +1 − (cid:1) possible choices. Letting N denote thelargest index for which ˜ k N (cid:54) = 0, the total number of edge arrangements is (cid:18) ˜ km (cid:19) = N (cid:89) n =0 (cid:18) ˜ k n m n (cid:19) . (3.5)Combining (3.4) and (3.5) yields a total tree count of { p ∈ S M | ( κ ( p ) , β ( p )) = ( k, m ) } = (cid:18) km (cid:19)(cid:18) ˜ km (cid:19) , as was to be shown. Theorem 3.1 Let ( τ (cid:48) , R ) correspond to an M -layer medium for some inte-ger M ≥ , and write | τ (cid:48) | = τ + · · · + τ M +1 . Then H ( τ (cid:48) ,R ) ( t ) = (cid:88) k ∈{ }× Z M + b ( R, k ) δ (cid:0) t − | τ (cid:48) | − (cid:104) k, τ (cid:105) (cid:1) where for each k ∈ { } × Z M + , the amplitude b ( R, k ) is given by the formula b ( R, k ) = (cid:88) ≤ m ≤ min { k, ˜ k } (cid:18) km (cid:19)(cid:18) ˜ km (cid:19) ( − R ) ˜ k − m R k − m T m + . (3.6) Proof . Given τ (cid:48) , observe that the arrival time for a transmission scatteringsequence p ∈ S (cid:48) M having κ (cid:48) ( p ) = k is precisely | τ (cid:48) | + (cid:104) k, τ (cid:105) . 17 May 2013 eter C. Gibson The combinatorics of scattering T occurring in the formula for the transmissionamplitude b ( R, k ) has the form T = (cid:113) − R (cid:113) − R · · · (cid:113) − R M , which is not a polynomial in the reflection coefficients ( R , . . . , R M ). So itturns out that b ( R, k ) /T is a polynomial in the reflection coefficients, while b ( R, k ) itself is not. (Recall that by contrast reflection amplitudes a ( R, k )are polynomial functions of the reflection coefficients.) (a) (cid:45) (cid:45) (b) (cid:45) (cid:45) (c) (cid:45) (cid:45) (d) (cid:45) (cid:45) Figure 8: The exact reflection Green’s function (a) for the example (4.1),and the underlying reflection coefficients plotted against two-way travel time(b). The exact transmission Green’s function (c) is plotted on a smalleramplitude scale, with the same reflection coefficients plotted against one-way travel time (d).The formulas presented in Theorems 2.1 and 3.1 make it easy to computethe reflection and transmission Green’s functions exactly up to a finite cutofftime T > 0. We illustrate this in the present situation by working out a17 May 2013 eter C. Gibson The combinatorics of scattering τ, R ) representing an M = 10 layermedium, with the following values. n Two-way travel timefrom z − to z n τ n R n . . − . . . − . . . − . . . . . . − . . . . . . . . . . . . − . . . . . . − . τ M +1 = 0, which corresponds to measuring the transmission re-sponse exactly at z M , formulas for G ( τ,R ) and H ( τ (cid:48) ,R ) given in Theorems 2.1and 3.1 were coded up in Mathematica. The resulting functions were com-puted up to a cutoff time of 5.38 seconds for reflection and 3.69 secondsfor transmission. (The respective computations took 3.37 seconds and 17.35seconds on a 1.8 GHz i7 processor.) The computed pulse trains are depictedin Figure 8. The reflected pulse train has 19 242 terms, and the transmittedone has 35 059 terms; the majority of these are of an amplitude which is toosmall to be visible in the plots.This example illustrates that Theorems 2.1 and 3.1, while derived bycombinatorial methods not usually associated with the analysis of PDEs,offer a computationally tractable representation of reflected and transmittedwaves. This provides a new tool for application to the various imagingmodalities where one-dimensional models are important. A Appendix: The underlying PDE Let u ( t, z ) denote the velocity (in the z -direction) of a material particle atdepth z and time t , and let p ( t, z ) denote the pressure. The medium evolves17 May 2013 eter C. Gibson The combinatorics of scattering ρ ∂u∂t + ∂p∂z = 0 (A.1a)1 K ∂p∂t + ∂u∂z = 0 . (A.1b)where ρ ( z ) is the material density at depth z and K ( z ) is the bulk modulus.For the sake of definiteness we focus on the velocity field u ( t, z ), although theresults can just as easily be formulated in terms of p ( t, z ). The initial con-ditions corresponding to a plane wave unit impulse propagating downwardfrom z − are u (0 , z ) = δ ( z − z − ) p (0 , z ) = (cid:112) K ( z − ) ρ ( z − ) δ ( z − z − ) . (A.2)Letting G ( t ) denote the (velocity) impulse response at z − , we have that G ( t ) = u ( t, z − ) ( t > , (A.3)the solution at depth z − to the system (A.1a,A.1b,A.2). Similarly, theimpulse response at z M +1 is H ( t ) = u ( t, z M +1 ) ( t > , (A.4)the solution at depth z M +1 to the system (A.1a,A.1b,A.2).The following facts are derived in [4, Chapter 3], [3, Section 2] andelsewhere. For 1 ≤ n ≤ M , let τ n denote the two-way travel time (for atraveling wave) across the n th layer of the above M -layer medium, and let τ denote the two-way travel time from depth z − to z . For 0 ≤ n ≤ M , let R n denote the reflection coefficient at depth z n relative to a wave travelingtoward the interface from above. Letting K n and ρ n denote the densityand bulk modulus inside the n th layer—with K − , ρ − and K M +1 , ρ M +1 denoting the respective values at z − and any point z M +1 below z M —thetravel times and reflection coefficients are given by the formulas τ n = 2( z n − z n − ) (cid:112) K n /ρ n and R n = √ K n ρ n − √ K n +1 ρ n +1 √ K n ρ n + √ K n +1 ρ n +1 , (A.5)for 0 ≤ n ≤ M . Note that − < R n < eter C. Gibson The combinatorics of scattering References [1] N. Bleistein, J. K. Cohen, and J. W. Stockwell, Jr. Mathematics ofmultidimensional seismic imaging, migration, and inversion , volume 13of Interdisciplinary Applied Mathematics . Springer-Verlag, New York,2001. Geophysics and Planetary Sciences.[2] L. M. Brekhovskikh and O. A. Godin. Acoustics of Layered Media I ,volume 5 of Springer Series on Wave Phenomena . Springer, Heidelberg,1990.[3] K. P. Bube and R. Burridge. The one-dimensional inverse problem ofreflection seismology. SIAM Rev. , 25(4):497–559, 1983.[4] J.-P. Fouque, J. Garnier, G. Papanicolaou, and K. Sølna. Wave prop-agation and time reversal in randomly layered media , volume 56 of Stochastic Modelling and Applied Probability . Springer, New York, 2007.[5] K. A. Innanen. Born series forward modelling of seismic primary andmultiple reflections: an inverse scattering shortcut. Geophysical JournalInternational , 177(3):1197–1204, 2009.[6] G. Kunetz. Quelques exemples d’analyse d’enregistrements sismiques. Geophysical Prospecting , 11(4):409–422, 1963.[7] J.-F. Le Gall. Random trees and applications. Probab. Surv. , 2:245–311,2005.[8] R. G. Newton. Inversion of reflection data for layered media: a review ofexact methods. Geophysical Journal of the Royal Astronomical Society ,65(1):191–215, 1981.[9] F. Santosa and W. W. Symes. Reconstruction of blocky impedanceprofiles from normal-incidence reflection seismograms which are band-limited and miscalibrated. Wave Motion , 10(3):209–230, 1988.[10] R. P. Stanley. Enumerative combinatorics. Vol. 2 , volume 62 of Cam-bridge Studies in Advanced Mathematics . Cambridge University Press,Cambridge, 1999. With a foreword by Gian-Carlo Rota and appendix1 by Sergey Fomin.[11] A. B. Weglein, F. V. Ara´ujo, P. M. Carvalho, R. H. Stolt, K. H. Matson,R. T. Coates, D. Corrigan, D. J. Foster, S. A. Shaw, and H. Zhang.17 May 2013 eter C. Gibson The combinatorics of scattering