A caustic terminating at an inflection point

We present an asymptotic and numerical study of the evolution of an incoming wavefield which has a caustic close to a curve with an inflection point. Our results reveal the emergence of a wavefield which resembles that of a shadow boundary but has a maximum amplitude along the tangent at the inflection point.


Introduction
In this brief paper we present an asymptotic and numerical study of a twodimensional wavefield with a caustic near the cubic curve y + 1  3 γx 3 = 0 as x → −∞ in the far field.The problem was introduced in [5] along with an explicit integral representation in [5, (4.7)] for the wavefield in a small neighbourhood of the origin, where the cubic has an inflection point and the caustic terminates.Unfortunately, [5, (4.7)] contained some errors, but also gave little insight into the outgoing field as x → +∞.This paper both corrects the errors in [5, (4.7)] and provides a detailed asymptotic description and a numerical solution of the outgoing wavefield.
The main rationale for this paper is that improved understanding of this wavefield may give intuition about the far-field solution of the famous Popov inflection point problem, which describes the scattering of an incident whispering gallery wavefield propagating along a concave portion of a scatterer boundary when it reaches an inflection point of the boundary.This canonical problem, which remains unsolved in closed form, has a lengthy history, the most recent review of which is [7].The Popov problem seeks a wavefield in the region y > − 1  3 γx 3 satisfying either a Dirichlet or Neumann boundary condition on the curve y + 1  3 γx 3 = 0, and tending to an Airy function close to this curve as x → −∞.In this paper we consider the evolution of an Airy function wavefield past the inflection point of the cubic in the absence of any boundary.
In §2 we introduce appropriate local curvilinear coordinates (S, N ) based on the cubic curve, and then derive a representation for the wavefield close to the inflection point in the form of an integral involving a complicated phase function of S and N , correcting [5, (4.7)].We then consider the outgoing field in the limit as S → ∞, which we show to have a transition between bright and dark in the vicinity of the positive x-axis.In §3 we reformulate the solution in terms of integrals studied in [4], which are appropriate for the numerical scheme introduced in [3], allowing us to present a validation of the asymptotic analysis.Finally, in §4 we offer some conclusions.

Asymptotic behaviour of the wavefield near the inflection point
We align a two-dimensional Cartesian frame so that, for sufficiently small values of x and y, the caustic has equation y + 1 3 γx 3 = 0, x < 0, for some positive constant γ.Local scalings on these coordinates, which both preserve this cubic profile exactly and allow us to model wave propagation through the inflection point at the origin are found to be x = k −1/5 X, y = k −3/5 Y , with X, Y = O(1).One can then seek an approximate solution to the Helmholtz equation (∆ + k 2 )ϕ = 0 of the form ϕ ∼ e ikx Ã (X, Y ) as k → ∞, with Ã satisfying the parabolic wave equation An alternative approach, adopted and reviewed in [5], is to switch to local curvilinear coordinates intrinsic to the curve y + 1  3 γx 3 = 0 involving arc-length s and normal distance n along and from it, respectively.It then follows that, correct up to and including terms of O k −1 as k → ∞, Thus, with the scalings ) and seeking an alternative local solution ϕ ∼ e iks A (S, N ), we have that, to lowest order, It then follows that (2.1) can be re-written in terms of S and N as which leads to the Popov equation In order to match with the incoming Airy function we require that where the first (−S) 1/3 term reflects the dependence of the incoming field on the curvature of the associated caustic as S → −∞ and the constants in the pre-factor multiplying the Airy function are inserted for later convenience.Since the data is analytic, we expect the solution of (2.2) to be analytic for all S, N .
In order to deal with the branch point at λ = −γS 2 , it is convenient to write This decomposition can be justified by deforming the integration path from one along ℑλ > 0 and showing that the contribution from the region near λ = −γS 2 is negligible as ℑλ ↓ 0. The remainder of this section will be based on (2.5).

The limit S → −∞
Noting that when S is large and negative, λ + γS 2 5/2 = −γ 5/2 S 5 1 + λ γS 2 5/2 , we can expand E for large negative S in I + .This gives that E ∼ 2N λ − λ 3 6γS + . . .and a formal asymptotic expansion shows that the leading term in which, by putting λ = (−4γS) 1/3 τ , is asymptotic to as S → −∞.Thus, since I − tends to zero as S → −∞, A matches with (2.3) as expected.An estimate of the rate at which the Airy function is attained for large values of −S would involve a complicated analysis which we will not attempt in this paper.

The limit S → +∞
The situation is more interesting when S → +∞, when it is convenient to write λ + γS 2 = −S 2 T for T ≥ 0 in I − so that dT. (2.6) where Φ = γN S 2 − γ 2 S 5 10 and is a measure of transverse distance from the x-axis.Similarly, writing λ+γS 2 = S 2 T for T ≥ 0 in I + leads to dT. (2.7) We begin by considering I − which we see from (2.6) has a stationary phase point at T = K.The contribution from this point to the value of the integral is exponentially small as long as S 2 K ≫ 1 and, in this case, the main contribution to I − comes from a region near T = 0 where T 5/2 is negligible to lowest order.However the T 2 term in the exponent needs to be retained to ensure convergence of the integral.Hence the leading term in the integral is where T = S 5 KT and ϵ = 1 S 5 K 2 .We now note that The device of introducing the parameter δ is one way of deriving higher order terms in asymptotic expansions of stationary phase integrals, as described in [6].It gives us the result that and it can also be used to obtain higher order corrections for the estimates we will derive in the rest of this paper.
The result corresponding to (2.8) for I + can be read off by replacing K by −K and it tells us that None of these results applies when K is of O(S −2 ) or smaller.In this regime, we write K = S −2 K and T = S −2 t, so that dt. (2.9) The dominant contribution as S → +∞ now comes from the stationary-phase point at t = K and we soon find I − is given by Moreover, the use of the regularisation leading to (2.8) reveals that the next order term in the stationary phase expansion of I − is of relative order S −1 .However, there is also a contribution from the end point t = 0, which can be estimated as in (2.8) to give a leading order term e −iΦ iS K .This can be shown to cancel with a corresponding term in I + in which there is no stationary phase contribution.Hence we conclude that, with an expected relative error of O(S −1 ), (2.11) as S → +∞ with K ≥ 0 and K = O(1).We note that the derivation of (2.10) from (2.9) involves the assumption that the stationary phase point K is sufficiently far from the origin that K ≫ O(S −1/2 ).We will examine the region where K = O(S −1/2 ) shortly.
The calculations above can be repeated to show that the contributions to both I + and I − near T = 0 cancel at least to leading order as S → +∞.Hence we need only consider I + , which has a stationary phase contribution.
It is now more convenient to write λ + γS 2 = S 2 τ 2 so that where (2.13) Hence there is always one stationary phase point in τ > 0, τ = τ 0 (K) say, which is the positive root of and is such that τ 0 (0) = 0 and τ 0 = O(|K| 1/3 ) as K → −∞.The leading order term in the expansion of the integral in (2.12) is thus , S → +∞, K < 0, (2.16) and so, as , which is equal to (2.11) when K = 0.However, the comments made after (2.11) apply equally to the derivation of (2.15) which only holds if S 5/2 τ 0 ≫ 1, and we will describe the inner layer when K = O(S −5/2 ) in the next subsection.
Meanwhile we note that the results (2.11) and (2.16) indicate that the farfield as S → +∞ is dominated by the region in which K is negative and of O(1).Plots of the numerical evaluation of A will be given in the next section and comparison will be made with (2.11) and (2.16).

K
In this region, the sign of K is no longer important and we will simply consider (2.7) and (2.6) when K = S −5/2 K. Writing T = S −5/2 t then, to leading order as S → +∞, (2.17) where v = K + t.Repeating the exercise for I − gives and hence Thus the amplitude |A| in this region is independent of K and matches with (2.11) as K → +∞ and with (2.16) as K → −∞.
Even though the amplitude of the far-field solution can be described analytically, we have only worked to the lowest order when obtaining the asymptotic expansions for A and this has resulted in A having discontinuous slope as a function of K. We expect the solution of (2.2) to be analytic everywhere and that the asymptotic approximation will become increasingly smooth when taken to higher orders but, as mentioned after (2.8), this is a challenging task.Hence it is helpful to compare these predictions with numerical calculations, which we do in the next section.

Numerical validation
In this section we compare the asymptotic approximations obtained in the previous section with accurate numerical evaluations of the integral (2.4), obtained using the PathFinder software [2].This implements the algorithm described in [3], which automates the numerical steepest descent method for oscillatory integrals (see e.g.[1, §5]), automatically performing appropriate contour deformations and dealing robustly with multiple coalescing stationary points.
To obtain an integral in a form amenable to evaluation using PathFinder, which in its current form applies only to integrals with polynomial phase, we apply a change of variable to rewrite (2.4) as A(S, N ) = −2e i(γN S 2 −γ 2 S 5 /10) Ã32 (S, N − γS 3 /3), where, adapting the notation of [4], with Γ 32 being any contour starting at t = e i9π/10 ∞ and ending at t = i∞.The tilde on Ã32 is included to indicate that the integral A 32 of [4] has been modified to include a non-trivial amplitude function F (t) = t.We note that Ã32 (X, Y ) solves the parabolic wave equation (2.1).
In Figure 1 we show a plot of |A| as a function of the curvilinear coordinates S and N , and in Figure 2 we show the corresponding approximate Helmholtz equation solution e iks(x,y) A(S(x, y), N (x, y)) as a function of the Cartesian coordinates (x, y).In Figure 3 we plot |A| as a function of N on the line S = −10, alongside the Airy function approximation (2.3).In Figure 4 we plot |A| as a function of K on the lines S = 5 and S = 10, accompanied by the approxi-  )).This reveals how, as S increases, the wavefield evolves from a beam-like structure to become more like a shadow boundary.In all comparisons we see excellent agreement between the asymptotics and numerics.

Conclusions
Motivated by the Popov problem of finding the outgoing wavefield generated by a whispering gallery wave as it approaches an inflection point, this paper addresses the wavefield that emerges when an incoming Airy function wave is centred on a curve whose curvature decreases to zero.An exact solution can be found in terms of the complicated integral (2.4) whose integrand contains a branch point.An asymptotic calculation using the stationary phase method reveals that the radiated field is strongest near the tangent at the inflection point but, unlike a Gaussian beam, it has an asymmetric structure whose largest amplitude is attained near this tangent and eventually resembles a shadow boundary.We have also computed the integral using the method of [3] after making a change of variable so as to obtain a representation without branch points.This has yielded wave profiles that compare favourably with the asymptotic predictions.