ENGINEERING RESEARCH INSTITUTE UNIVERSITY OF MICHIGAN ANN ARBOR SOLUTION OF A COMBINED REFLECTION AND REFRACTION PROBLEM BY USE OF A TRANSPORT EQUATION W. C. MEECHAM Project M936 OFFICE OF NAVAL RESEARCH, U.S. NAVY DEPARTMENT CONTRACT N6 onr-23221 June 25, 1952

INTRODUCTION In the treatment of propagation of radiation in an inhomogeneous medium, it is frequently necessary to consider both reflection and refraction effects. A method of treating problems of this type where some of the radiation leaving the bounding surface is refracted back to it is set forth here. The method depends upon two major assumptions. First, it will be assumed that the components of radiation contributing at a particular point have random phase. Secondly, it will be assumed that ray-tracing techniques may be used for any particular component. A transport equation is then derived whose solution yields the measurable quantities. Solutions are obtained at large ranges which express the intensity at any point in terms of the velocity structure, The solution at large distances is essentially an equilibrium of the radiation field. ii

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN SOLUTION OF A COMB/INED REFLECTION AND REFRACTION PROBLEM BY USE OF A TRANSPORT EQUATION The basic problem to be considered is the following. We are given a half-space bounded by some irregular surface which does not deviate too far from a plane. The surface is free so that the boundary condition is that the velocity potential (0) in the case of an acoustic field vanishes at the surface. The velocity of propagation in the medium depends on depth (z) alone. We will further assume that the sources are localized. A special class of problems in which we will be particularly interested is one in which the velocity structure is such that some rays which start near the surface will be bent back and return to the surfaceo The region in which this occurs is frequently called the channel. Further let it be assumed that changes in the velocity structure are slight in a distance comparable to one wavelength of the radiat ion (x). Mathematically the problem would be stated0 find a function $(xy,z) which fulfills the following equation and boundary condition: [72 + k2(z] 0 = 0o where k2 is a function of z alone. 0 = 0, when z = p(x,y), where p is the function which describes the free bounding surface. Further, we are given some system of sources which gives rise to the radiation field. It is immediately evident that an exact solution of the above problem is beyond present methods. Even if we: consider the two simpler derived problems, one of pure reflection and one of pure refraction, only special problems in these subclasses lend themselves to anything like exact solutions. The introduction of the assumption that the function k2 is slowly varying will be the key assumption which leads to the method to be described herein. If the function k2 has this property we will be able to speak Qof rays and be able to trace their course through the medium by the familiar "ray tracing" techniques. As mentioned previously, we will

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN be particularly interested in the special class of problems in which some rays can return to the surface when they have been started near the surface. This is a situation roughly similar to optical focusing. z p(x) Fig. 1 We will now collect all the principal assumptions to be used in the construction of the transport theory (see Fig. 1). a) For the moment we will consider the two-dimensional problem in which all quantities depend on x and z alone. b) The incident radiation will consist of plane waves propagating in a direction such that their propagation vectors are perpendicular to the generating element of the (cylindrical) free surface p(x). c) Further let it be assumed in the original formulation that k2(z) is real. This implies that there is no attenuation in the medium. d) We will be particularly interested in the solution far enough away from the source of the plane waves so that those rays which are returned to the surface have "bounced" many times. e) We will assume that the phases of the radiation are random, so that we can calculate the total intensity of the radiation at any point by adding up the intensity of the individual contributing waves. In fact, in the formulation we will consider only the intensities of the contributing components.

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN f) Finally we assume that k2(z) is slowly varying in a length x. Let us pause here and briefly review ray-tracing theory, deriving those functions which will be used later, Suppose we' wish to solve the wave equations ( V'2 + k2)0 = O under the approximation that k2 is slowly varying in a length k. Let us then substitute 0 = X(xy,z)'e-'iWW(xYZ)/Co dropping terms under the approximation stated above, we arrive at1 (-V )2 = C2/c2, which is often called the eik;onal equation. Here c is the velocity at a point x,y,z, and co is *the velocity at a reference point. Now if we further restrict our problem so that k2 is a function of z alone (the velocity a function of z alone), we can easily deduce the following relationship: cos e c cos 0 co (See Fig. 2.) Fig. 2 Fig. 2 1. See,. for instance, Physics of Sound in the Sea, Part 1: Transmission, National Research Council (abridged).

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN We can now find the range, x, of a single ray making an angle g0 with the:.surface, as a function of the depth of the ray, z, from the following equation, making use of the stated relationship for cos @: x =f cot c dy = cosso c c dy x = o or 0 = Once the function c is given we can find the distance which a ray starting from the surface must travel before returning to the surface. This distance, given for all initial angles, constitutes the function L(Qo) of this work. Let us now define Z(%o,z) as the ratio of the intensity of a ray bundle at a depth - z to its intensity at the surface; this ratio is, of course, different for ray bundles moving at different angles %0. In order to find this function, imagine- a plane wave starting from the origin and confinedbetween %0 and 0o + dOo. Suppose the source is radiating F units of energy per radian per second in the direction 0o Now, since no energy is lost from the ray bundle, its intensity at P will be given by I - Fi o ds where ds is the perpendicular distance between:the two confining rays. It is clear from the diagram that ds = -.. sin @ do; hence we have F - 1/ - sin 0. F =o Now after using the equation for the range, c cy x = cos 9J A '. a We arrive at the following expression Z e= l/c02 sin 0o sin c c dy " '. -0 — L: _ ':,,: A:

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN Let us now simplify this expression, Suppose that c/co f(z) Then let z = g(c/c0). Now make a change of variable in the above integral, substituting c = c/co. We arrive at Z = l/sin g0 sin / ) ( If we now integrate by parts, we have C, g,,~, f(z) Z =./sn1./s in cos 'o(1 2 os2 ) 1/2 f(z) ~ d_~g (1 o2 cos2 o)1/2 cos2 2OJ This integral can be carried oit nunerically, if desired. One case of interest is that for which the velocity is a linear function of dept;h, z.o ln this Case the integral in. the denominator vanishes. Let us now return to the derivation of the transport equation. Before writing down the general equation describing the system, let us define the.functions to be used and state -the viewpoint which leads to the formulation. First let us conseider what we will call a statte "",. By this we will mean a bundle of rays which makes an angle 0 with -the positive x direction when leaving the surface (see Figo 1). It will be assumed, as already mentioned, that the velocity structure is such t:hat we can use ray-tracing techniqules. The entire effect of the velocity structure is contained in a function L(O) which is the distance between consecu.tive bounces of a state 90 L(9) is symmetric about 0. =/2 (.ee Fig0 3). i I |gn, is the maximum value of 0 which I I stays in the charnel. 0 GM T-M Fgo, 3 The effect of the rough surface will be placed in a fu.nction. A.,( 9 ). The:function A gives the intensity of the plane wave reflected into

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN the direction 9 due to a plane wave of intensity one incident upon the surface from direction @'., It has two properties which Will be useful in the development. a) A(@, G' ) = A(',G)., This property follows from the reciprocity theorem. It is completely general and holds for any surface. b) A(@,@')d%. = 1. This merely expresses the fact that no energy leaves the medium upon reflection at the bounding surface. We will consider a state G at a distance x from the origin 0 placed at the source. (Something should be said here about the source of the radiation, In general, of course, the source will not give rise to a single plane wave. We will analyze the radiation from any source or system of sources as a sum of plane waves~ ) We will add up the amount of intensity lost from this state 0 by collisions with the surface and the contributions to this state from collisions of other states with the surface. Call the intensity of the state @ at x, I(Gx)o I measures the total intensity of states G at x. It must be remembered that x and G do not completely describe the radiation. We might say that there still remains a degeneracy. At a given position x there can be a great number of states, all of which belong to. (see Fig. 4). of the states shown only the group A will collide with the surface at the position x. However, the sum of the intensities of all of these groups (their intensities will be assumed equal, which is consistent with the previous assumptions) makes up the functioan I(g,x)o Thus, in order to find the intensity of the states 9 which collide with the surface at a given position between x and x+dx, one must take I(',x)ax L(G) From.these considerations we are' led to the following equation governing the function I:.. = (O A(GXv) d( A(',))d' (1) L(G").(.) 6~~~~~~~~~~

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN z A C D Fig. 4 The first integral. is the contribution to 9 from all of the states Go colliding with the surface and being scattered in the direction.. The second integral measures the loss of intensity from the state 0 due to scattering from the surface~ The factor 1/L(Q), as explained above, is needed to give: the probability that a given state 9 collides with the surface in a given range interval. Now usinrg the property Equation (1) reduces to: + I(9. l I(',x) A(,')d' (2 x L() L() (2) This integro-differential equation is the one which governs the function I. Let us now discuss a method of solution of Equation (2) and some general properties of its solution. We will use here a separation of variables technique in order to obtain a solution of (2). With this in mind, let I(G,x) = e()X(x), where as usual G is a function of 9 alone and X is a function of x alone. If we assume that the e~ form a complete set we arrive at the following solution to (2): I(ex) =2ccze-x )C(e)

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN where O() satisfies:. cc((.) +.() - L(',) A(G,G')d@'. (3) L(G) L(') This process reduces the integro-differential equation to a homogeneous integral equation. We will now continue with a discussion of some properties of the solutions of (2) and (3). First let us prove that the ce's involved in Equation (3) are all real. To do this return for the moment to (2)2. Consider the Fourier expansions of all quantities. Since the angular range is (O,-), we will need to consider only even terms in the expansions. Let, I(x, 9) Zdn(x)e2inG l/L () ne2inG = ()e2inG n — - Then, substituting in (2): ' n(x)e + 2in n(x)e2inG n= m -; n=m:P m(-te )e2mG k n(x)e2inG d.O 2in@4nx 2imG mn=- 0 2 For a similar discussion, see M. C. Wang and E.o Guth, Physo Revo, 84, 1092 (1951)o

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN or DX (x)e 2nG + K Z r-n n(x)eir r= -Vrn=where D = - dx From this it follows, because of the linear independence of the functions e2in that: D (x) +) L / r(X) Inr( )d rs r o or r0 O We have here an infinite simultaneous system of homogeneous linear differential equations. Now let us take a, "core" frm this system, for instance by limiting the summation of r between + Ro We have then, R I Fnr D +dnr qJnQnr (')d'?QLJr 0= r=-R O Now in this case we know that each tVr(x) is a solution of the symbolic equationl i D + i-j ij ('91 )d 1 = O. 0 Here| | means the determinant of the matrix. This is a linear homogeneous differential equation with constant coefficients (in general of order 2R + 1). The solutions are exponential, of type e-~x, where the allowed c.g's are roots ofo || Cij D + ti-j jj i )dB | Oy (4) with D considered an ordinary variable, But the roots of this equation are just the negatives of the e igenvalues of the matrix:

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN M.-ij =jijj - ~i-j(Gf)dQ We will show that this matrix is Hermitian, in which case the roots (W. 8s) are all real. Consider the matrix Mo If n = n* and.n = q n* the matrix is Hermitiano However, L(G) and A(9,GO) are realX Therefore their expansion coefficients have the above property. We thus conclude that all the permissible acs in Equation (3) are realo One is also tempted to conclude from this treatment that the allowed. acs take on discrete values, no matter how many equations from the full infinite system were used. Of course this is not. a proof. We will now prove that the) appearing in Equation (3) are orthogonal on the interval (O0,.) with the weighting function l/L(9).o From Equation (3) we see: _m1 o (o) + ( Ore() a::(.)l (.,) 0)' -'2~O () + Oc(G) = _ _2(G) A(o,G)dG' L(G) L() Now multiply the first by (e3()/L(G), the second by ( )/ (o integrate with respect to 9, and. subtract0 The resulting equ tion iso 0(a sL(G) d.' - I Oca2 ( )O.zl) Cl ( ) A(9, ) dGd.B' o 0- "()L(o') (5) /O 0 O ( 9)()A(O) gde, o o. (e)L(e' ) However, since.A(O9',G) Awe see that the right side of (5) is zero, Thus o dG = O if Cl C2' (6) ~L(9) 10

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN which is the required property. We will now show that B93 -6; 0 a(G)dG 0 if o O O 0 To do this, integrate Equation (3) with respect to 0. The result is: -.:8 + G,()do = e o:.~' ) o L(g) OA(0 0 L(G' o )Ad.'dG Now if we can interchange integrations on the right side (it is assumed throughout.that the functions are sufficiently regular to permit such operations ) we have cO0 = 0, (7) where the normality condition on A has been used. This was the required relationship, 'We will now show a property which I(0,x) must,possesso (Actually this property follows from Equation (7) above ). Consider (2) again Jx + L ( 9 ) J, A (.:G@' ) dQO' ~'(O~ ' L() A(,G )dG ' Integrating with respect to @ leads to,x.) J. + 0 (G) -,:( A(@Gx)d )gd + O ~x od:r..,..(Q"): r.., (.' ), 'd' Using the normality relation on AT(e,,),, we arrive at o I(@e,x)dG = Io Thus this integral is constant in range. This may at first seem like a contradiction, since some of the beam is being lost frcm the region near the surface as we progress in range o However, it must be.remembered that I(gx) measures the total intensity of the state @ at the position x, and it is immaterial whether or not the ray happens to lie near the surface at this poaition. 11

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN We will now outline a plausibility argument for the orthogonality and reality proofs. Equation (3) can be written in the form: o L('i) (@) ( t d@ Suppose we change the independent variable. We arrive at _@a(0 my - _ X(0) where. = - 1/ca; de. = d4/L.('), or - = 0(9); and the bar over the function means it is expressed in terms of 0 instead of 0. This equation could be written: OL(0) = X Q (0')K(0,n0)d.i a,b,K being defined in the obvious way. Note that K(0',"') is symmetric. Now if K were continuous on the interval [b], we could deduce the following from the Hilbert-Schmidt theory: a) X's (aCs) are real, b) AIs are discrete, c) OG: are real (except perhaps for multiplying constants), d) 0Qc are orthogonal, and e) any function which could be written as:f(x) = K(x,t)g(t)dt (8) a where g(x) is continuous on [ab], could be expressed as a linear combination of the eigenfunctions of the kernel K. Now it must be remembered that K is not continuous on [ab We might, however, hope that a sufficiently good approximation to the solution could be obtained by taking an approximate delta function which was continuous on the range. Itn this case all the above would follow. Even in this case, we could not hope to solve every boundary condition 12

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN which might be given but only those which could be expressed in the form shown in Equation (8). If it should turn out in any particular problem that the eigenfunctions are complete, there is no difficulty. If not we can do one of two things. First we might hope to approximate the initial radiation distribution (this is the boundary condition) with the eigenfunctions at hand. If this is not good enough we would have to use one of the more elaborate schemes which have been developed for this type of problem. 3 Let us now consider the solution of the physical problem in more detail. We will assume that the G form a complete set. We had I(oGx) = 4Sce- ca(O). (9) The coefficients cac are determined by the boundary conditions, For instance, suppose the distribution in 0..at x =.0 is p (g)o Then we have p() =:ab(3 cc..from. which A~ 1 1. (p,.(),d C - _ ___ L(@) (@)do after using (6). This relationship is sufficient to solve the problem. Frequently we will not need the complete series. For instance; 1. If we are interested in the region x = 0 to +oo, then all negative values of CX must be excluded. 2. If further we wish the result at large ranges (as is often the case) only the smallest values of CX are important. In fact, at great ranges the important term has Ca = 00 This assumes that the al's are discrete. 3. In the above situation, when x>> L(G) 3 MC. C. Wang and E, Guth, loc, cit. 13

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN I(x,G) = co 0(~). Now (o(0) may depend on A(G,G') and L(G), but not on x. It determines the intensity variation with depth. It seems, then, that the intensity pattern at various depths does not change at different ranges in this approximation. Let us adapt the solution to experimental conditions. We must consider two complicating effects. A) In the typical problem we would have a localized source located in the channel. In order to correct for the geometric (cylindrical) spreading we should multiply the result by l/x. B) If k2(,z) is slightly complex, there will be a slight damping. This is usually the case experimentally. To correct for this, introduce an exponential factor e-Px where X is the parameter describing the attenuation. (1/p is, of course, the distance at which the intensity of a plane wave is reduced to 1/ of its initial value.) The final result is: I(R8)- e-P R(0 cO(0), R where we have replaced x by R. From this and ray-tracing diagrams, we can immediately derive the intensity variation with depth. In order to do this we will need one more function which is readily obtainable by ray-tracing techniques (see above). This is a function which we have called Z(G,z). Z(G,z)/R is the intensity of a ray @ at depth z and range R. Also note the probability of finding a ray G at depth z, and range R is I(R,Q)/L(@). Then if }'is the intensity, =,d Z(Oz:,B (() r~= e-i~R Co T/ a(@) Notice that the integral is essentially the coefficient of the term ac = 0 in the expansion of Z(@,z) in the system of functions kU(.). co is, except for a constant, the total energy radiated. This can be seen by using the condition (7) or (9). For then, 14

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN / o.p(G)dG c - oo()dG It is not difficult to see from () that one solut for is L(@) Whenever L(G) becomes very large for some G, as is the case in the ordinary channel problem, the constant c becomes small. If now we return to our equation for A, 'R '-= e-R coc S z(GZ) dG, (10) Using the expression already developed for z we have here the expression for the intensity at large ranges in terms of the velocity-depth function. This, then, is the final result. There are many refinements in viewpoint and treatment which can be made. To name a few: 1. By utilizing solutions other than ca = 0O we could obtain predictions of the intensity closer to the source. 2. A more careful treatment of the function 0o(0) is required 3. It would be interesting to examine in greater detail the effect of directivity and depth of the source on the intensity at large ranges. 4 The question of uniqueness is an interesting one. Under certain conditions it is possible to get other solutions. In order that (3) have a unique solution when O =., it is necessary that A(GGO) be non-reduced in the group theoretical sense. In other words, for uniqueness it should not be possible by permuting rows and columns to put A in a form consisting only of smaller matrices located along the diagonal (all other elements being zero). 15