## Abstract

A saddle point formula for integrals whose integrands possess a multi-valued exponent, and one or more poles is derived. The result includes a number of error (or Fresnel) functions, which is equal to the product of the number of poles and the number of sheets possessed by the exponent's Riemann surface. Such formulae have previously been suggested on the basis of known exact results and transformations, both of which are specific to a particular exponent. In general, the multi-valuedness can be taken into account by a straightforward modification to a standard procedure, which would ordinarily yield only a single error function for each pole. In the context of diffraction theory, this type of approximation remains valid in the vicinity of an optical boundary, even in some cases where a wavefield is incident in a direction almost or exactly parallel to a sharp obstacle (grazing incidence).

## 1. Introduction

The purpose of this article is the derivation of an improved saddle point formula for integrals whose integrands possess poles and a multi-valued exponent. For ease of presentation, and physical interpretation, we consider a class of such integrals, members of which occur frequently in diffraction theory, for example in the half plane problem illustrated in figure 1. Thus, let(1.1)where the location of the pole *α*=*α*_{0} is given by(1.2)with *k*>0 and *Θ*∈[0,*π*). The exponent *Χ* is defined as(1.3)with(1.4)so that *ϕ* is a solution to the Helmholtz equation in polar coordinates (*r*, *θ*), i.e.Integral solutions to other partial differential equations, which have different exponents, can be treated in a similar manner. The contour consists of the real line, traversed from left to right, with indentations above *α*=−*k* and *α*=*α*_{0} and below *α*=*k*, see figure 2, below, and the branch of the exponent is chosen so that *γ*(0)=−i*k*. This corresponds to the situation in which *ϕ* is derived from a time-dependent wavefield, with harmonic factor e^{−iωt}. The integral in equation (1.1) generally represents the outgoing diffracted field generated when the plane wave term strikes one or more sharp obstacles. We are concerned with the evaluation, in the far field *r*≫1, of the contribution from the first order saddle point *α*=*α*_{s}, where(1.5)and, in particular the effect of the pole *α*=*α*_{0} upon this procedure. The function *f* is taken to be free from essential singularities; this is generally the case for integrals arising in diffraction problems. In addition, for the purposes of the derivation we shall assume that it has no poles. We generalize our result to deal with multiple poles, and those outside the interval (−*k*, *k*) at the end of §2. If *f* should possess branch points other than *α*=±*k*, then these may contribute to the solution in some cases. We therefore choose to present our discussion in terms of the method of steepest descents, rather than that of stationary phase, since if the path of integration has to be diverted around a branch cut, one can often show that the resulting contribution is exponentially small. In addition, certain standard integrals appear automatically when the calculation is carried out using this technique.

One way to view the method of steepest descents, which is particularly effective in problems involving poles, is a contour deformation, followed by a complex mapping, which transforms equation (1.1) into a real line integral to which Watson's lemma may be applied (for large *r*), provided that the saddle is not close to the pole. This results in a non-uniform approximation, whose greatest merit is simplicity; it forms the basis of many results in ray theory (Keller 1962), however, it becomes singular as the saddle approaches the pole. Alternatively, following the transformation, the pole can be removed by the addition and subtraction of a simple integral which can be evaluated exactly. The result is a uniform approximation, consisting of a scaled complex error function (Abramowitz & Stegun 1965; henceforth, we shall just write error function for brevity), and an integral which is approximated by Watson's lemma. The former is *not* an approximation; it is a component of the exact value of *ϕ*, and satisfies the governing partial differential equation. It accounts for the rapid, but continuous, activation and deactivation of the residue term as the variation in *θ* causes the saddle to cross the pole *α*_{0}. In place of the error function, many authors use the complex Fresnel function (Noble 1988), however the formulae that follow are slightly more compact when given in terms of the former. The relationship between the two is given in §2. The non-uniform approximation can be retrieved by replacing the error function with the first term in its asymptotic expansion for large argument. The uniform approximation has been applied to the solutions of diffraction problems involving Helmholtz equations by several authors, including Rawlins (1975) and Abrahams & Wickham (1990). For a general discussion of these methods, the reader is referred to Bleistein & Handelsman (1986) and Felsen & Marcuwitz (1994). As we shall see, even the uniform approximation breaks down in certain limits. The main result of this article is the fact that the multi-valued nature of the exponent function *Χ* in equation (1.1) causes the steepest descent mapping to generate as many poles as the Riemann surface possesses sheets (two in this case). We derive an improved uniform asymptotic formula for (1.1) which consists of two error functions and, as usual, an integral to be approximated by Watson's lemma. In order to give a physical interpretation to these remarks, we now consider the integral (1.1) in the context of the diffraction problem shown in figure 1, in which a wave is incident at angle *Θ* upon a barrier occupying the half line *θ*=0. This may represent diffraction of acoustic waves by a half plane (Noble 1988), or of bending waves by a crack in a thin plate (Norris & Wang 1994) for example. The resulting scattered field can be derived via the Wiener–Hopf technique and expressed as an integral of the form (1.1) (additional terms, representing evanescent modes are also present in the latter case). Typically, the function *f* in (1.1) takes different forms on either side of the half line *θ*=*π*, though the field is of course continuous here. The merger of the saddle and pole occurs at the optical boundaries of shadow and reflection shown in figure 1. This ought to be expected, as these divide the regions in which the incident and reflected fields are present from those where they are not. Typically, two separate uniform approximations must be derived, one that includes an error function to account for the activation and deactivation of the incident field, and one that does likewise for the reflected wave. The former is applied for *θ*≤*π*, and the latter for *θ*≥*π*. (In a problem with a more complicated geometry, it may be necessary to consider several different regions.) However, in the limits *Θ*→0 and *Θ*→*π*, which are known as ‘grazing incidence’, the optical boundaries are close together, and in their vicinity, both uniform approximations break down. It is reasonable to suppose that, in fact, both error functions are present in the solution for all *θ*. However, to directly establish this, we would need to assume differentiability to all orders, which is certainly not apparent across *θ*=*π* given the integral form of *ϕ* (equation (1.1)). Various alternative approaches exist for integrals satisfying the Helmholtz equation. One such method is to transform the integral via the so-called ‘complex angle substitution’(1.6)This results in a steepest descents problem in which all functions are single valued, though the exact evaluation of contour integrals in the *β* plane is not straightforward (Noble 1988). Another is to deduce the improved approximation by comparison with the exact evaluations of the Sommerfeld half plane integrals, as in Borovikov & Kinber (1994). Finally, were we to make the assumption that the field is indeed differentiable to all orders across *θ*=*π*, then a ray method could be employed, as in Lewis & Boersma (1969). For integrals satisfying other differential equations, which have a different exponent function *Χ*, however, there may be neither exact results in literature, nor a convenient mapping such as (1.6). The present method can be applied in these cases, even when the exponent function takes a complicated form, as it does in problems involving wave propagation and diffraction in anisotropic media, see for example Norris & Achenbach (1983), Felsen & Marcuwitz (1994) and Thompson & Abrahams (2005). In addition, it explicitly shows the precise means by which multiple error functions arise from the asymptotics of diffraction integrals possessing only one pole.

## 2. Derivation of the improved approximation

Return to equation (1.1), consider *θ*∈(*θ*,*π*), and deform the contour onto the path of steepest descent . This interval for *θ* is sufficient to illustrate the method; the approximation in the limiting cases *θ*=*0* and *θ*=*π* will be determined from our final results by continuity. Also, we must have(2.1)here, since the residue from *α*=*α*_{0} eliminates the incident plane wave term in the shadow region *θ*<*Θ*. The expression for *ϕ* now becomes(2.2)where *H*(.) represents Heaviside's unit function. In the case where *α*_{s}=*α*_{0} (i.e. *θ*=*Θ*), is indented to the right, leaving the pole beneath the contour and hence we take *H*(0)=0. Given that(2.3)it is not difficult to show that the steepest descent path is defined by the equation(2.4)where *u*≥0, and the different roots represent the two branches emanating from the saddle point. The arcs at infinity yield no contributions, and we can use equation (2.4) to show that the descent path always crosses the real line precisely twice, at the saddle and also at the point *α*=−*k* sec *θ*. The branch points *α*=±*k*, therefore cause no difficulty. It should be emphasized that equation (2.2) is continuous across the half line *θ*=*Θ*; the effect of the contour crossing the pole accounts for the presence of the Heaviside function. A typical steepest descent path in the *α* plane is shown in figure 2. The angle at which the contour crosses the interval (−*k*, *k*) turns out to be of some importance. To determine it, we first note that(2.5)where the prime symbol denotes differentiation with respect to *α*; this convention is maintained throughout. It is not difficult to see that this expression is always positive imaginary at the saddle point, in fact(2.6)If we now writethen in the limit *α*→*α*_{s}, we have(2.7)note that the left-hand side of this expression is clearly positive real for . Given that , and therefore , moves from upper to lower half plane as (−*k*, *k*) is crossed, it follows that the tangent to the descent path at the saddle point intersects the real line at a clockwise angle 3*π*/4 (see figure 2). Let us now introduce the steepest descent mapping(2.8)with(2.9)In the *t* plane, the (first order) saddle is fixed at the origin, and the image of the steepest descent path is clearly the real line, since the exponential in the integral is now simply . The location of poles in the *t* plane, the direction in which the path of integration is to be traversed and the manner in which it should be indented if necessary must all be carefully considered. The function *g*(*t*) cannot usually be determined for general *t*, however its value (or perhaps residue) at certain critical points, such as the origin, can be found relatively easily. The location of poles in the *t* plane can be deduced from the following observations. First, we note that the integrand in equation (1.1) exists not on a single complex plane, but on a Riemann surface, which possesses at least two sheets, due to the branch points in the exponent. If the function *f* also possesses branch points (as is generally the case), then there are often multiple sheets. On each sheet of this Riemann surface, there is a simple pole located at the point *α*=*α*_{0}. Since d*t*/d*α* is finite, it follows from equation (2.9) that the function *g*(*t*) has poles at corresponding points in the *t* plane. Those poles which lie on sheets with *γ*(0)=−i*k* are mapped by equation (2.8) to the point *t*=*t*_{0}, which is defined up to a factor ±1 by(2.10)whereas those on sheets with *γ*(0)=i*k* are taken to the point *t*=*t*_{1}, with(2.11)The next step is to determine the direction in which the real line is to be traversed in the *t* plane, and any appropriate indentations that may be required. To achieve this, we first note that, except on indentations, for . Second, since we must take a square root to obtain *t*, we avoid the point *t*^{2}=0. Thus, in the *α* plane, is indented to the right of the saddle point in a semi-circular arc with radius *ϵ*. This is consistent with the indentation used previously in the special case *θ*=*Θ*; it now applies for all *θ*. The radius *ϵ* must be sufficiently small that the orientation of with respect to the pole *α*=*α*_{0} is not changed. On the arc, we may apply equation (2.7), with *η* decreasing from an initial value of 3*π*/4 to −(*π*/4) (see figure 2*a*), to show that arg(*t*^{2}) varies from 2*π* to 0 as the indentation is traversed. It now follows that the image of in the *t*^{2} plane traverses the lower side of from right to left, encircles the origin, and finally traverses the upper side from left to right, as shown in figure 2. If we now take the square root , we find that the real line is traversed from left to right in the *t* plane, and that the pole *t*=*t*_{1} lies on the line arg(*t*)=*π*/4, i.e. it is always above the path of integration, henceAs usual, the surd symbol is used only when taking the positive root of a positive, real quantity. The situation with regard to *t*_{0} is slightly more complicated, since for *θ*=*Θ*, we have *t*_{0}=0, i.e. as *θ* varies, this point may ‘jump’ from one sheet of the square root's Riemann surface to the other. Consequently, we cannot appeal to any complex analysis argument to determine its position, however this is unnecessary—continuity of the resulting expression in *θ* will suffice. We now let the indentation radius *ϵ*→0 in all cases except *θ*=*Θ*. Clearly, *t*_{0} lies on one side of the contour for *Θ*>*θ* and on the opposite side otherwise; only this behaviour can maintain continuity given the presence of the Heaviside function in equation (2.2). The real line has an indentation above the origin in the critical case *θ*=*Θ*, therefore we may conclude that *t*_{0} lies on the upper side of the contour for *Θ*>*θ*, and on use of equation (2.10) we obtainWriting(2.12)and separating into partial fractions, we find that *ϕ* is now reduced to the form(2.13)where(2.14)(2.15)and(2.16)In the first instance, the path of integration *Γ* passes below the pole in all cases, and a residue has been included where appropriate. Since there is now no possibility of a pole crossing an integration path, the critical value *h*(*t*_{0}) follows immediately; we must havein order to maintain continuity in *θ*. This result can also be obtained by using equation (2.12) in (2.9) and rearranging to yieldthe limit *t*→*t*_{0}, which implies *α*→*α*_{0} may now be taken without difficulty. The procedure described above cannot easily be used to determine *h*(*t*_{1}), since in this case the function *f* must be evaluated on a sheet which does not contain . If *f* has branch points (as is usually the case), then it is not always obvious which sheet should be chosen. Thus, at this stage, we simply writewhere *c* is a constant to be determined. In fact, we shall see later that *c* is the reflection coefficient for the barrier in question. Finally, we determine the value of *h*(0); thus we differentiate equation (2.8) twice with respect to *α*, to obtain(2.17)Next, on substitution of equation (2.12) into (2.9) and evaluation at the point *t*=0, which corresponds to *α*=*α*_{s} we find that(2.18)where we have used the identity *α*_{s}−*α*_{0}=i*t*_{0}*t*_{1}. The argument is determined from equation (2.17) by noting the direction of variation of *α* and *t* as the paths of integration pass through the saddle in the respective complex planes. As *t* passes from left to right through the origin *α* travels down a path whose tangent intersects the real line at a clockwise angle 3*π*/4 (see figure 2*a*), and hence . It now remains to simply collect together all of our results, evaluate the standard integrals (2.14) and (2.15) and take the leading order behaviour of (2.16) via Watson's lemma. It should be emphasized that this is the *only* approximation applied at any stage of the entire procedure. After some manipulation, we find that(2.19)(2.20)and(2.21)Here, we have made use of the standard integral (Abramowitz & Stegun 1965)in fact this is valid for all *z*_{0} provided that the path of integration is adjusted so as to pass below the pole. The notation *w*(.) represents the scaled complex error function, i.e.this is related to the aforementioned complex Fresnel function *F*(.) viaAs noted above, the terms *ϕ*_{1} and *ϕ*_{2} are not approximations; they are exact solutions to the Helmholtz equation with the special property that they may ‘contain’ a plane wave. To see this, we need the two results (Abramowitz & Stegun 1965)(2.22)and(2.23)Thus, *ϕ*_{1} includes the incident field for *θ*>*Θ*, however the plane wave included in *ϕ*_{2} cannot be activated for *θ*∈(0, *π*). This is as we should expect; there are no such terms other than the incident field in this region. Indeed, the only other plane wave in the entire problem is the reflected field, which is present for *θ*>2*π*−*Θ*, precisely the region in which the plane wave in *ϕ*_{2} is activated. In fact, from equation (2.22), we can writeand immediately we recognize the plane wave term as the reflected field. The constant *c* is thereby determined as the reflection coefficient for the barrier occupying *θ*=0. More rigourously, *c* can be obtained by repeating the entire procedure for *θ*∈(*π*, 2*π*) and imposing continuity as *θ*→*π*. The two simpler approximations can be retrieved by using equations (2.22) and (2.23) in (2.13). In particular, the non-uniform approximation, with residues omitted, is given by the first term in equation (2.21), i.e.(2.24)The effect of the additional terms is to regularize this expression at the two poles *θ*=*Θ* and *θ*=2*π*−*Θ*, and to include the plane waves where appropriate. If we now consider the integral (1.1) and the approximation given by (2.13) with (2.19)–(2.21) as functions of the parameter *Θ*, an immediate generalization to deal with integrals possessing poles outside the interval (−1,1) becomes available. Thus, for a simple pole located at *α*=*α*_{p}, we write *α*_{p}=−*k* cos *Θ*_{p}, where . Addition of the correction terms and , where(2.25)to (2.24) now yields the improved formula. The coefficients *A*_{+} and *A*_{−} can be determined by comparison of residues—the resulting approximation must be regular at *θ*=*Θ*_{p} and *θ*=2*π*−*Θ*_{p}. In both cases, we can determine whether to take the upper or lower sign by ensuring that the first term in (2.22) is activated in the correct region (cf. equation (2.23)). Thus, we have reached the result suggested by Borovikov & Kinber (1994), via a method which does not depend upon comparisons to other diffraction integrals whose values are known a priori, and can therefore be applied in cases, where the governing equation is not that of Helmholtz. The procedure of correcting the diffraction coefficient can simply be repeated to account for integrals whose integrands possess multiple poles. For *Θ*_{p}∈[0, *π*], equation (2.22) yields plane waves as above, whereas if , we obtain evanescent modes (these may propagate without loss in some special direction, such as along a barrier, and can therefore play an important role in the far field behaviour).

Before examining a more complicated problem, it is useful to briefly investigate the properties of the approximation given by equations (2.13) and (2.19)–(2.21) by applying it to the classical problem of diffraction by a half plane. The geometry is identical to that shown in figure 1, and the boundary condition requires that ∂*ϕ*/∂*θ*=0 for *θ*=0 and *θ*=2*π*. The solution to this problem in the form (1.1) is (Noble 1988)(2.26)where the fractional power in the denominator is defined so that for *α*>−*k*. The reflection coefficient *c* is easily shown to be unity, and on recalling that equation (2.13) is derived for *θ*∈(0, *π*), the third term now vanishes identically. We are left with(2.27)which clearly does not break down in any limits—it is the exact solution to the problem in question, originally obtained by Sommerfeld (1896) (for an English reference, see Sommerfeld (1954)). Thus, in this case nothing is actually *gained*; we have simply retrieved an old and widely known result. It is worth noting, however, that had we taken *h*(*t*_{1})=0, we would obtain a formula in place of equation (2.27) which is singular in the limit *θ*+*Θ*→0. Another point of interest is that the improved formula continues to give the correct result for all *θ* if we multiply (2.26) by csc *Θ*/2 and then set *Θ*=0, thereby merging the pole with the branch point of the function *f*.

## 3. An example from thin plate theory

Finally, we consider the situation in which *ϕ* represents the transverse displacement of an infinite thin plate, with a rigid strip along the half line *θ*=0. Then *ϕ* satisfies the partial differential equation (Timoshenko 1940)with boundary conditions for *θ*=0 and *θ*=2*π*. The solution to this problem, originally derived by Norris & Wang (1994), is(3.1)where we have introduced the conventions that a function shown without argument depends upon *α* alone, and a subscript ‘0’ denotes evaluation at the point *α*=*α*_{0}. It should be emphasized that our notation is rather different from that used by Norris & Wang (1994). Here, *γ* is defined by equation (1.4), and the second exponent by(3.2)with for . The contour is the same as before; it passes between the branch points *α*=±i*k*. A superscript ‘+’ (‘−’) denotes a function that is analytic on and above (below) ; these are obtained via an analytic product factorization (Noble 1988) so that *K*^{+}*K*^{−}=*K*, etc. The function *K*(*α*) is given byand *K*^{+} can be expressed in terms of finite, non-singular contour integrals (Norris & Wang 1994). Thompson & Abrahams (2005) apply a procedure that requires only one such integral to a similar problem. The remaining factorizations can be performed by inspection, thusandwhere the branches are chosen so that , and . It is important, for the purposes of our subsequent analysis, to briefly discuss the possibility of evanescent modes in equation (3.1). First, consider terms with an exponent involving *λ*. A straightforward calculation shows that these have a saddle point located at *α*=−i*k* cos *θ*, and that the associated steepest descent path is defined bywhere *u*≥0. Clearly, this is the real line if cos *θ*=0; it lies in the upper or lower half plane for cos *θ*<0 and cos *θ*>0, respectively. In the latter case, a diversion around *α*=−*k* is necessary, note that the terms in question have no branch point at *α*=*k*. All contributions are exponentially small, except in cases, where *θ*≈0 and *θ*≈2*π*, i.e. in the vicinity of the rigid strip. Next, consider the terms with an exponent involving *γ*. These have a branch point at *α*=−i*k*, but not at *α*=i*k*. We can apply the method of steepest descents as before; equation (2.4) shows that the deformed integration path crosses the imaginary axis at the point . Therefore, it must be diverted if cos *θ*>|sin *θ*|, or else it will pass below *α*=−i*k*. Since *γ* is pure imaginary for *α*=(1+*u*)i*k*, *u*>0, the largest contribution due to the diversion comes from the branch point itself, and it decays exponentially, at least as rapidly as . To summarize, there are evanescent modes present for all *θ*, though these do not contribute significantly to the far field in general. There are also modes that propagate along the faces of the rigid strip *θ*=*0* and *θ*=2*π*; these are evanescent in other directions. If we now treat equation (3.1) as an integral of the class (1.1) (i.e. disregard terms with exponent *λ*), with(3.3)we should expect a good approximation, except possibly in the immediate vicinity of the rigid strip. The reflection coefficient may be obtained by taking the residue in (1.1) at *α*=*α*_{0}, on the lit side of the strip; we obtain(3.4)Introducing a subscript ‘s’ to refer to evaluation at the saddle point, we find that, for *θ*∈(0, *π*),(3.5)it is generally simplest to perform manipulations in terms of *γ*_{0}, *λ*_{s}, etc., where possible. Equations (2.13), (3.4) and (3.5) define an approximation that is valid for all *θ*; the same result occurs if we repeat the entire procedure considering instead *θ*∈(*π*, 2*π*). For *θ*=0 and *θ*=2*π*, equation (3.5) vanishes identically, due to the factor . Interestingly, this means that, if *Θ* is sufficiently large to permit the use of (2.23), then (2.13) accurately predicts the displacement at the strip (zero). In addition, if *Θ*=0 (head-on incidence), then *c*=−1, and equation (2.13) satisfies the boundary condition exactly. Thus, the omitted modes, all of which are evanescent except possibly along the axis of the strip, make a significant contribution to the displacement here for small, non-zero *Θ* only. Finally, it is of some interest to assess the accuracy of equation (2.13) for *θ*=*π*, since the uniform approximation is known to break down in the vicinity of this half line as *Θ*→*π*. Thus, return to (3.1), set *θ*=*π* and deform the path of integration to pass along the faces of the branch cuts in the upper half plane. Ignoring exponentially small contributions, we find that(3.6)wherein(3.7)and(3.8)Here, the first integral has no pole at the point *α*=*α*_{0}, and the subscript ‘R’ has been introduced to denote the fact that *γ*^{−} is to be evaluated on the right face of the branch cut. Following an application of Watson's lemma in equation (3.7), the substitution reduces both integrals to standard forms, which evaluate to giveSetting *θ*=*π* in equation (2.13), and making use of the identities (2.2) and *c*−1=2*λ*_{0}*K*_{0}, we can derive precisely the same expression for *ϕ*, showing that the improved approximation is valid here for all angles of incidence.

## 4. Concluding remarks

We have demonstrated the derivation of an improved saddle point formula for integrals possessing poles and a multi-valued exponent. For ease of presentation, the procedure is performed upon an integral which satisfies the Helmholtz equation; it can easily be adapted to other cases. The resulting formula includes a number of error functions, which is equal to the number of poles multiplied by the number of sheets possessed by the exponent's Riemann surface. The method is easy to apply, regardless of the number and location of the poles, and the range of validity often extends beyond the interval for which it was derived.

## Footnotes

- Received February 16, 2005.
- Accepted November 25, 2005.

- © 2006 The Royal Society