Logarithmic catastrophes and Stokes's phenomenon in waves at horizons

Waves propagating near an event horizon display interesting features including logarithmic phase singularities and caustics. We consider an acoustic horizon in a flowing Bose-Einstein condensate where the elementary excitations obey the Bogoliubov dispersion relation. In the hamiltonian ray theory the solutions undergo a broken pitchfork bifurcation near the horizon and one might therefore expect the associated wave structure to be given by a Pearcey function, this being the universal wave function that dresses catastrophes with two control parameters. However, the wave function is in fact an Airy-type function supplemented by a logarithmic phase term, a novel type of wave catastrophe. Similar wave functions arise in aeroacoustic flows from jet engines and also gravitational horizons if dispersion which violates Lorentz symmetry in the UV is included. The approach we take differs from previous authors in that we analyze the behaviour of the integral representation of the wave function using exponential coordinates. This allows for a different treatment of the branches that gives rise to an analysis based purely on saddlepoint expansions, which resolve the multiple real and complex waves that interact at the horizon and its companion caustic. We find that the horizon is a physical manifestation of a Stokes surface, marking the place where a wave is born, and that the horizon and the caustic do not in general coincide: the finite spatial region between them delineates a broadened horizon.


Dedication
This paper is dedicated to Sir Michael Berry in celebration of his 80th birthday. Of his many contributions to physics and physical asymptotics, one of the major themes of his work is his interest in waves near singularities, ranging from the most dramatic occurrences such as tsunamis [1,2] and tidal bores [3,4], to the most gentle (yet profound) in the form of Stokes's phenomenon [5,6,7,8]. In particular, he has devised minimal models for undular bores that reveal in a characteristically clear way the central role played by caustics as well as an analogy to the Hawking effect [3], and it is a related connection we pick up here in the context of a flowing superfluid. Starting in the 1970s [9] and continuing today [10], Michael Berry has championed the application of catastrophe theory to caustics and we humbly follow in his footsteps in this paper, focusing on novel catastrophes with logarithmic phase singularities that accompany an acoustic event horizon in a superfluid.

Introduction
If the flow speed of a fluid exceeds the speed of waves in the fluid then the latter are unable to propagate against the flow and an analogue of an event horizon occurs. This situation is capable of mimicking many aspects of black hole physics where the effective spacetime metric depends on the fluid flow and acoustic waves play the role of light [11,12,13,14,15]. Analogue event horizons for classical waves, including the classical analogue of Hawking radiation, have been observed in water tank experiments [16,17,18,19,20,21,22] and also in optical fibres [23,24,25]. Another system where analogue event horizons can be created is a flowing Bose-Einstein condensate (BEC) formed from ultracold atoms. These are among the simplest examples of superfluids and are so cold that quantum processes such as the analogue of spontaneous Hawking radiation can dominate thermally activated phonons. Analogue Hawking radiation in BECs has been anticipated theoretically for over twenty years [14,26,27,28,29,30,31,32,33,34,35], and was recently realized in a series of experiments by J. Steinhauer and coworkers [36,37,38,39,40,41]. In the present paper we focus on event horizons in BECs where excitations obey the Bogoliubov dispersion relation, but many of our results apply to event horizons in general.
Tidal bores (shock waves that travel up rivers due to a rising tide being funnelled into an V-shaped estuary) are a dramatic wave phenomenon that also involves hydrodynamic event horizons, as pointed out by Michael Berry in his studies of undular bores [2,3]. His work emphasizes that the wave front corresponds to a caustic where two waves coalesce, and this is the approach we adopt in this paper. The connection between event horizons and caustics is illuminating because caustics take on certain universal shapes described by catastrophe theory that are structurally stable against perturbations (and hence generic), and furthermore each catastrophe is dressed by a unique wave function. In the simplest case of two coalescing waves the caustic is a fold catastrophe and the wave function is the Airy function [42]. References [2,3] show that the wave profile of an undular bore can indeed be expressed in terms of an integral which has the Airy function as its kernel.
In this paper we find that in the vicinity of a horizon in a BEC the waves satisfy a third order differential equation in position space which is similar to the one obeyed by undular bores. However, we generalize the treatment given in [3] to the case where the wave frequency ω is nonzero and show that this leads to an Airy-like wave function but with an additional logarithmic contribution ω log k to the phase, giving a "log-Airy function". A logarithmic term has previously been shown to arise in aeroacoustic flows from jet engines [43], in certain path integrals in radio astronomy as studied by Feldbrugge et al from a wave catastrophe point of view (albeit one which differs from both our own and the standard approaches) [44,45], and also for gravitational black holes when dispersion is present [46]. The fact that a logarithmic phase singularity is a very general feature of waves in accelerated frames (such as near event horizons) has been highlighted by U. Leonhardt and collaborators [47,48], who first suggested a connection between horizons and wave catastrophes in a general sense, although a characterization in terms of the Thom-Arnold catastrophe theory which underlies standard caustics was not pursued. Links between caustics, Airy functions, and black hole physics have also been previously reported in more detail by Nardin et al in the context of water waves with subluminal dispersion [49]. They pointed out the presence of a fold catastrophe bifurcation by taking a 'dynamical systems' approach which is similar in spirit to ours. However, we will show in this paper that the logarithmic phase term connects the wave structure near the horizon to higher wave catastrophes than the fold.
Indeed, the main motivation of the present paper is the observation that the coefficient multiplying the logarithmic term in the log-Airy function adds a second control parameter and thus the underlying catastrophe has some of the character of a cusp catastrophe, the next in Thom's hierarchy of catastrophes beyond the fold (which only has a single control parameter). This feature is supported by classical solutions of Hamilton's equations describing the motion of wavepackets which show that an event horizon gives rise to a broken pitchfork bifurcation in (z, k) phase space where a single solution bifurcates into three. The universal wave function dressing a cusp catastrophe is the Pearcey function, but unlike the Pearcey function the logarithm leads to mathematical complications in the evaluation of the integral representation of the solution, such as the need to introduce branch cuts and multiple Riemann sheets. We describe a systematic method for handling these complications using exponential coordinates [43]. Along the way we shall identify instances of Stokes's phenomenon which occurs when an exponentially small wave appears behind a dominant wave and constitutes the "quietly beating heart of asymptotics" [50]. In particular, we find that the horizon itself is the edge of a Stokes surface, giving some physical meaning here to this mathematical concept.

Bogoliubov dispersion relation
Elementary excitations in a BEC obey the Bogoliubov dispersion relation [51,52] where c = m ∂n ∂P −1/2 is the speed of sound expressed in terms of the mass m of the atoms and the compressibility of the gas ∂n/∂P (n is the number density and P the pressure). k c = mc/ is analogous to the Compton wavenumber (inverse of the healing z 0 supersonic subsonic Figure 1: Schematic representation of a quasi-one dimensional acoustic black hole formed within a flowing fluid. The blue circles represent the constituent particles of the fluid, and the blue arrows indicate their direction and magnitude of flow. The vertical dashed line at z = 0 indicates the position of an acoustic event horizon so that the flow is supersonic when z < 0 (inside the black hole), and subsonic when z > 0 (outside). The wavepackets represent Hawking radiation (pairs of back-to-back phonons) that propagate away from the horizon. One phonon escapes and moves rightwards, while the other gets swept leftwards into the black hole.
length) [14,28], and provides the characteristic scale at which the Bogoliubov relation becomes dispersive: at small wavenumbers such that k k c the relation reduces to the linear form ω = ck but for k > k c the relation curves upwards and hence is supersonic, or "superluminal" in the gravitational context. This implies that there are no true event horizons in BECs because there can be waves with arbitrarily large speed. However, these are energetically suppressed and, crucially, there is still a bifurcation when the flow speed exceeds the speed of sound so that essentially the same phenomena arise as in fluids with subluminal dispersion relations [16,49,53,54]. In a Lorentz invariant system the dispersion is purely linear, however this leads to divergences at the horizon known as the trans-Planckian problem which can be resolved by the introduction of dispersion [12,14,28,46,55,56,57,58,59,60,61,62,63,64,65].
A schematic representation of a localized event horizon in a fluid is shown in figure  1. An early suggestion for experimentally realizing such a situation was to consider a BEC flowing around a ring that has a localized constriction such that the flow speed inside the constriction is forced to increase above the speed of sound in order to maintain the current and there is no buildup of atoms [26]. This set-up actually generates two horizons: a black hole horizon where the fluid enters the constriction and from which long wavelength phonons cannot escape and a white hole horizon on the other side where phonons cannot enter. In the eventual experiments performed by the Steinhauer group a potential step (formed by a laser) is swept through a long thin cigar-shaped BEC and this effectively creates a waterfall over which the atoms flow [36]. We shall not attempt to model the details of these situations but instead content ourselves with the simplest theoretical model and consider a quasi-one dimensional BEC with a stationary velocity profile u(z) that close to the horizon changes linearly in space [53,63] u(z) ≈ −c + κz (2) where the velocity gradient κ (analogus to the surface gravity of a graviational black hole) is positive and c is the speed of sound at the horizon. This describes a flow from right to left whose magnitude is linearly increasing in the direction of the flow. At the point z = 0 the flow speed is the same as the speed of sound and hence this is the location of the horizon: for z > 0 the flow is subluminal and for z < 0 it is superluminal (the speed of sound in a BEC made of a dilute gas of ultracold atoms is of the order of millimetres per second so this is easy to achieve).
The frequency that appears in Eq. (1) is that in the fluid's rest frame, but an observer in the laboratory frame (where the horizon is fixed in space) will see a Doppler shifted frequency. A Galilean shift between the two frames gives the relation ω = ω − u(z)k where here and from now on we use the unprimed notation to represent the laboratory frame. Thus, ω will be the frequency (energy) of the excitations seen from the laboratory frame. It is constant if the flow is time independent and must obey ω ≥ 0 in order to be physical. Conversely, ω can be negative and this leads to an analogy to antiparticles.

Nondispersive case: trans-Planckian singularity
We will first consider the long wavelength nondispersive regime where the nonlinear terms in the Bogoliubov dispersion relation can be ignored. The linearized dispersion relation in the laboratory frame is given by In the limit where the flow varies very little over the wavelength of the produced sound waves, one can use a ray picture (geometrical acoustics). Treating ω(z, k) in Eq. (3) as a hamiltonian function, Hamilton's equations read These rays describe the centre of mass motion z (t) and k (t) of wavepackets [14,28]. Eqns. (4) can be numerically integrated in time starting from an initial choice of position and wavenumber (z i , k i ), and in the left image of figure 2 we have plotted the trajectories for two different initial conditions, one just inside and one just outside the horizon, the idea being that these show the fate of spontaneous excitations starting near z = 0. We see that in the phase space provided by the canonical variables (z, k) the horizon behaves like a hyperbolic fixed point: rays move away from the horizon and tend asymptotically to z = ±∞ whilst the magnitude of their wavenumbers are reduced asymptotically to zero (red shifting). These two solutions can be interpreted as the classical manifestation of Hawking pair production where one excitation escapes the black hole and the other is trapped inside [53]. Indeed, both solutions have the same (positive) value of the lab frame frequency ω but have different signs of the fluid rest frame frequency ω = ω−u(z)k and hence can be considered to be a particle-antiparticle pair. If we instead run the time evolution backwards we find that these solutions of the linear model track back to the horizon where they develop wavenumbers of infinitely large magnitude ±k. This is the trans-Planckian problem where the horizon seems to connect low energy physics to infinitely high energy physics. with a finite value of ω > 0. The plot is in the (z, k) phase space describing the position and wavenumber of wavepackets. The BEC is flowing from right to left: the right hand side of this figure is outside the black hole, and the left hand side is inside it. The red and blue branches φ out ±k are classical analogues of Hawking pairs that propagate away from the horizon towards z = ±∞. The wavenumbers diverge at the horizon. Right: Same as the leftmost plot but now for a nonlinear dispersion [Eq. (8)]. The resulting broken pitchfork bifurcation is the structure associated with a cusp catastrophe. Both solutions start at the initial point z i , and the critical point z c marks the turning point of the blue branch and is the location of a caustic where two real solutions merge to become two complex rays. Of these, only the decaying complex ray φ out ↓ is physical and is denoted by the gray dashed line.
We have only included the positive roots of Eq. (3) in figure 2. The negative root also gives a physical solution, however, it describes waves that move with the flow which pass through the horizon relatively undisturbed and so will not be included here, not least because they are almost completely decoupled from the Hawking effect [46,63]. However, including these other solutions is important if one wants total momentum i k i and total energy i ω −u(z)k i conservation in the fluid rest frame (the laboratory frame frequency ω is automatically conserved by Hamilton's equations).
In the left plot of figure 2 we have introduced the notation φ out k and φ out −k to label the two solutions. This is adapted from the paper by Coutant, Parentani, and Finazzi (CPF) [63] upon which much of the present work is based. In particular, we use a wave scattering formalism where "out" indicates that a wave is moving away from the horizon (in either direction), whereas the subscript ±k indicates whether the wave has a positive or negative wavenumber. In the dispersive problem we treat in the next section we will also have "in" waves that propagate towards the horizon. Finally, we note that within the linear theory the two sides of the horizon are disconnected: trajectories cannot cross the horizon.

Wave theory in the nondispersive case: logarithmic phase singularity
To obtain the effective wave equation that describes these waves we will follow Berry [3] and treat the frequency ω(z, k) as a hamiltonian which can be canonically quantized: ω(z, k) →ω(ẑ,k) where z →ẑ and k →k = −i∂ z . Using Eq. (3) we find that Schrödinger's equationωψ = ωψ for this system takes the form (the question of operator ordering in this equation is addressed in [3,63]). This equation describes a stationary solution ψ(z, ω, t) = ψ(z, ω)e −iωt where and A is a constant. The logarithmic phase singularity at the horizon means the phase is undefined there and is a further manifestation of the trans-Planckian problem. Note that we have explicitly included the eigenvalue ω within the argument of ψ to later make connections to catastrophe theory.

Classical solutions in the dispersive case: broken pitchfork bifurcation
In terms of the laboratory frequency ω, the full Bogoliubov dispersion relation is We will again keep only the positive root and furthermore will expand the right hand side and keep only the first nonlinear term. This is sufficient to capture the essential physics (a caustic where two or three waves coalesce depending on whether ω = 0 or ω = 0, respectively). We therefore work with the cubic dispersion relation The trajectories this generates via Hamilton's equations are plotted in the right image of figure 2. The initial conditions (z i , k i ) are indicated by black dots, and this time we choose the initial position z i for both trajectories to be inside the horizon as this allows a fuller picture of the dynamics: if the positive k branch (red) were to start outside of the horizon it would continue to propagate to the right away from the horizon, whereas here we see that it passes through the horizon. Thus, the Bogoliubov dispersion relation allows the connection of the inside and outside of the black hole. In fact, comparison of both the left and right plots in figure 2 shows that the effect of nonlinearity is to generate a broken pitchfork bifurcation such that for some values of z inside the horizon there are three different possible values of k. The structure of the solutions means the divergence of the wavenumbers at the horizon is now eliminated but despite this change the Hawking process is robust to the introduction of dispersion because at large length scales we still have a pair of outgoing solutions φ out ±k and only the near-horizon behaviour is strongly modified [12,46,55,56,57,58,59,60,61,62,63,64,65].
Like in the linear case, the negative k solutions (blue) in the left plot of figure 2 correspond to antiparticles because ω < 0 when k < 0. Their dynamics restricts them to be inside the horizon but this time they can travel in either direction. We recall that the motion of wavepackets is dictated by the group velocity dω/dk = dω /dk + u(z) and bearing in mind that in our setup u(z) < 0 this means dω/dk can be positive or negative irrespective of the sign of k, although for a large enough magnitude of k it will always be positive due to the superluminality of the Bogoliubov dispersion. The critical point z c shown in the right plot of figure 2 marks the place where φ in −k turns around and becomes φ out −k which is a caustic because turning points are places where two solutions coalesce. This coalescence "interacts" strongly with the third solution provided by the positive k branch to give the broken pitchfork structure. This is the signature of a cusp-type catastrophe which has two control parameters: one determines where we are in the bifurcation (the z coordinate) and the other determines the degree to which the pitchfork is broken (this will turn out to be ω). Right at the cusp point itself, which occurs at ω = 0 and z = 0 (see figure 4), we find an unbroken pitchfork where all three solutions coalesce. The caustic at z c also marks the point where two complex rays appear. One is physically reasonable since it corresponds to a decaying mode φ out ↓ . The other corresponds to a growing mode φ out ↑ so we do not include it in our considerations [63].

Length scale associated with quantum effects
The wave function for the dispersionless case given in Eq. (6) is scale free and has no intrinsic length scale to cut off the logarithmic divergence. However, the presence of k c in the dispersive case allows the introduction of a new length scale which provides the characteristic distance over which the horizon becomes smeared out [46]. For the remainder of this work we will work in units where lengths are scaled by d, and times are scaled by κ −1 : The length d has a quantum origin because it depends on the Compton wavelength h/mc which sets the scale below which pair creation becomes important. The role that d plays in regulating the logarithmic phase singularity suggests that it is 2nd quantization of the classical wave theory that ultimately resolves our analogue trans-Planckian problem even though we do not explicitly 2nd quantize in this paper. In the classical theory of the pair production, the nonlinearity is sufficient to regulate it.
3.6. Wave theory in the dispersive case: log-Airy function Canonically quantizing Eq. (8) (in scaled units) yields a third-order linear ordinary differential equation in z-space This approximately describes ψ(z, ω) for both small z and small k, although we keep in mind that although k has been approximated to be small, it is not so small as to be equivalent to the linear dispersion approximation ω−u(z)k ≈ ck. Ignoring the nonlinear dispersive effects equates to removing the ∂ 2 z operator within the brackets, and Eq. (11) reduces to the linear case of Eq. (5) as expected.
Although Eq. (11) may be solved for ψ(z, ω) exactly in terms of a linear combination of three 1 F 2 hypergeometric functions, it is more instructive to use a Fourier transform to give an equivalent integral representation [46,63] where A is an arbitrary complex constant. Note that although Eq. (11) is only valid for small k, we can extend the domain to include all k ∈ (−∞, ∞) by assuming that the integrand in Eq. (12) primarily contributes near k = 0 and vanishes as k becomes larger in both the ±k directions. With careful choice of cuts, this turns out to be true in our case. The wave function in Eq. (12) is the primary object of interest in this work. We refer to it as the log-Airy function because the cubic polynomial in the exponential in the integrand is analogous to that of an Airy function [66,67] (see chapter 9 of [68] for the standard modern definition), which is the universal wave function that dresses a fold catastrophe [69,70,71,72], but differs due to the parametrically prefactored additional logarithm in the phase. There is a temptation to combine the logarithm with the pole term to produce a complex order branch cut, however, as we shall see, it proves analytically important to retain this term in the exponent.
Eq. (12) has appeared before in a number of related contexts. The most closely related to our present treatment is the study of black hole radiation in Lorentz-violating theories in the papers by Coutant and Parentani [46] and CPF [63], where the coefficient multiplying the logarithmic term is inversely proportional to the surface gravity which in turn sets the Hawking temperature. Similar logarithmic catastrophes have also been studied by Stone et al [43,73] in the context of aeroacoustic noise, although Eq. (12) differs from [43,73] in that it contains a 1/k pole term within the integrand, and its integration range is from k ∈ (−∞, ∞). An integral of the form of Eq. (12) also appears in Berry [3] in the context of tidal bores. Berry works in a frame of reference that moves along with the bore such that the bore is stationary, which is equivalent to setting ω = 0. This collapses Eq. (12) back to the integral studied by Boyd [74] without the logarithm within the phase.

Strategy
Our goal now is to investigate the behaviour of Eq. (12) by characterizing its behaviour in the (z, ω) plane through a study of saddlepoint contributions and Stokes's phenomenon that determine the physical ray and wave behaviour. The negative k range of the integration leads to differences (and simplifications) from the work of Stone et al [43,73]. In particular, the complex logarithm within the phase requires a choice of branch cut in order to define ln(k) when k < 0. Following CPF [46,54,63,75], there are two physically motivated choices of branch cut which can be made for our system in the complex k-plane: (i) along the negative imaginary k-axis, leading to ln(k) = ln |k| + iπ when k < 0, (ii) along the positive imaginary k-axis, leading to ln(k) = ln |k| − iπ when k < 0.
We will refer to (i) as the +iπ choice of branch cut, and (ii) as the −iπ choice. Through large |z| asymptotics of Eq. (12), CPF determine that not all the waves scattered by the horizon can be obtained by a single choice of branch cut, and a description involving contributions from both is required. In fact, a complete solution of the third order differential equation Eq. (11) requires three linearly independent global modes obtained by different choices of branch cuts and contours. Luckily one of the modes (the one associated with the exponentially growing wave φ out ↑ ) is not needed for the Hawking problem and we will ignore it here.
Our approach differs from CPF because we use exponential coordinates that allow us to focus on the near horizon behaviour in greater detail. In particular, CPF deform the contour of integration in the k-plane of Eq. (12) towards steepest descent paths over the saddles. However, some of the resulting deformed contours snag on the branch cuts and do not intersect saddles. In order to stay on the principal Riemann sheet, they make the approximation that ω is small and pull the logarithm within the phase out of their asymptotic approximation (see equation 57 in [63]). To deal with the snagged contours, they use an additional approximation referred to as the "dominated convergence theorem" which requires that k c → ∞ (so that the dispersion is approximated to be linear). These two separate approximations result in a partial loss of the near horizon wave structure, although it accurately describes the asymptotic behavior sufficiently far from the horizon. By contrast, exponential coordinates allow us to treat multiple Riemann sheets in a transparent manner and we freely pass through the branch cut onto other sheets where we pick up additional saddlepoints in a straightforward way, without the need for additional asymptotic approximations. Furthermore, our method allows us to resolve caustics (coalescence of saddlepoints) and Stokes's phenomenon (sudden change of contour as new saddlepoints appear) which are the elemental processes that determine the structure close to the horizon.
We shall treat the lab frame frequency of the emitted quanta ω as a real parameter. For the system in question ω is small and positive (to observe the Hawking effect). The details of the calculation will be explained for this case, although a similar approach can be taken for both larger ω and ω < 0.
4.2. The large parameter, self-similar scaling properties and the classical limit Saddlepoint methods work well in the semiclassical regime where the phase oscillates rapidly (apart from at the saddlepoint itself). A large parameter λ that performs the role of an inverse Planck's constant −1 can be introduced into the log-Airy function by changing variables in Eq. (12) (13) Like in the cases of the standard diffraction integrals (Airy, Pearcey, etc.), this procedure gives redefined control parameters, in this case z/λ 2/3 and ω/λ. Therefore, in addition to the redefinitions given in Eq. (10), we henceforth make the further redefinition The powers (2/3, 1) to which λ is raised are known as Berry indices and determine how the fringe spacing in the wave function evolves in the directions specified by the control parameters as λ is changed [76]. The index 2/3 matches that of the Airy function, but the unity power for ω does not have a counterpart in the standard diffraction integrals, see table 36.6.1 in [68] (the Pearcey function has the indices 3/4 and 1/2). Another difference to the standard diffraction integrals is that the scaling does not change the overall magnitude of ψ(z, ω) since the new factor e −i(ω/3) ln λ outside the integral is a pure phase term. This is different to the Airy case for which the amplitude diverges as λ 1/6 , where the exponent 1/6 is the Arnol'd singularity index [68,70]. For the Pearcey function it is 1/4. Thus, the classical limit λ → ∞ leads to an infinitely rapidly varying phase of the wave function but does not lead to an infinite amplitude as it does for standard caustics.
What physical parameter should we choose for λ? The natural choice is the quantum length scale d defined in Eq. (10). More precisely, we choose the dimensionless ratio where d 0 is an arbitrary reference length scale. The classical limit in our problem is therefore d 0 /d → ∞ where pair creation occurs only at vanishingly small length scales and the spectrum is linear to infinitely large values of the wavenumber. Recalling that z in Eq. (13) is already scaled by d, we find that the redefined coordinate in Eq. (14) grows as d −1/3 . In other words, the fringe spacing in physical space shrinks, as expected.
Our approach to the large parameter needed for the saddlepoint analysis is different to that employed in other works such as that by CPF where their parameter (equation 34 in [63]) is spatially dependent and vanishes at the horizon. By contrast, our λ is spatially constant which is important because we seek to resolve the saddlepoint structure even at the horizon.

ω = 0 case and pole contribution
To make a link with previous work, we first study the case when ω = 0, which is related to a zero-energy "soft mode". In that case the logarithmic term in Eq. (13) vanishes and it reverts to an integral which appeared in Berry [3] in the context of tidal bores. The integrand then becomes effectively an Airy-function exponent with a pole and may be evaluated as an integral of the Airy function itself and was studied by Boyd [74].
This integral contains three asymptotic contributions, two from the saddlepoints and one from the pole, see figure 3. There is but a single caustic point at z = 0, where the two saddles coalesce with the pole at s = 0. As z runs from z < 0 to z > 0, the steepest descent contours in the s-plane will deform and reconnect as is the case for the Airy function. As the system passes through z = 0 the contours unavoidably cross the pole and so generate an additional residue contribution, essential for the step function of the initial data in the case of the bore. The overall result is the well known bore form of a step function modulated by an Airy function.
By contrast, we shall see for ω = 0 below that although the waveform at fixed z mimics that of a bore, the rise in the overall magnitude for z < 0 is not due to the pole. Rather an equivalent analysis shows that steepest paths in the s-plane do not ever cross the pole at s = 0. They deform around it, but never generate a residue from it. The modulated step function appearance of the exact result can be seen to be generated from pairs of saddlepoint contributions, where one of the saddles comes from adjacent Riemann sheets.

Similarity to cusp catastrophe
We now study the case of ω = 0. The location of the saddlepoints on the principal sheet of the k-plane is given by ∂ s f (s, z, ω) = 0. For fixed ω = 0 the locations of the saddles s j are therefore given by Hence the presence of the logarithm increases the number of saddlepoints by one to three, rather than the two underpinning the Airy function. Thus, we expect the analytical skeleton of the log-Airy function Eq. (13) to be more akin to that of the next most complex function in the catastrophe theory hierarchy, the cusp, whose waveform is given by the Pearcey function [68,69,71,72] Pe(y, x) = Ψ Cusp (y, x; 1) .
The diffraction pattern generated by this integral is depicted in the left hand panel of figure 4 and features a two dimensional set of fringes generated by three wave interference below a cusp-shaped caustic (black curve). The highest intensity fringes are close to the caustic which is the place on which the saddles of Eq. (17) coalesce in pairs except right at the tip of the cusp where all three saddles coalesce. Asymptotic expansions about the saddlepoints diverge on the caustic but the exact waveform given by Eq. (17) is smooth: interference between the waves resolves the ray singularity on the caustic. The overall wave behavior of the Pearcey function relative to the caustic (depicted by the black curve) is shown in the right panel of figure 4. The red curve gives the position of the Stokes set for the Pearcey function (as first found by F. J. Wright [77]). The ray limit of the Pearcey function, which provides the scaffold upon which the waves are hung, is obtained from the saddlepoints of the exponent of the integrand in Eq. (17) [77], where an evanescent wave is born. Right: The number of waves (saddles) in each region relative to the cusp caustic and the Stokes set. The notation "wave" is short for "real wave" (i.e. non-evanescent).
log-Airy function with (x, y) playing the roles of (−ω, z). This implies that, at least as far as real rays are concerned, the ray structure for the log-Airy function is identical to that of the cusp. The position of the caustic for the log-Airy function is given by the condition for coalescence of two or more saddles by the simultaneous satisfaction of ∂ s f (s, z, ω) = 0 and ∂ 2 sf (s, z, ω) = 0. Eliminating s from these two equations yields Since the real saddles of the log-Airy function obey the same equation as those of the Pearcey function the classical caustic structure for the two functions is also identical. In fact, the local similarity to the next highest form in the catastrophe hierarchy of functions will be true for all such cuspoid diffraction integrals perturbed by logarithms in their exponent [68]. Despite the identical caustic structure for the log-Airy and Pearcey functions, their waveforms differ significantly. It is a relatively straightforward numerical calculation to evaluate the log-Airy integral in Eq. (13) and the results for λ = 1 for both choices of branch cut are displayed in figure 5. The two dimensional fringe pattern of the Pearcey function inside the cusp caustic and its damped nature outside are not replicated for the log-Airy function. Rather, for values of z below the caustic the log-Airy function has one dimensional fringes. For the same range of ω these are more pronounced for  Although we will only study these specific points, the respective regions they occupy also display equivalent behavior for ω > 0. For ω < 0 the dominance of contributing saddles is simply the opposite of the corresponding ω > 0 regions, but we only focus on ω > 0 in this paper.

Steepest descent contour diagrams
We perform a careful steepest descent analysis of Eq. (13), extending the work of CPF [63] to study particular regions of interest relative to the horizon and caustic. Due to the fact that the horizon and the caustic do not coincide except at ω = 0, there is effectively a broadened horizon (gap between horizon and caustic) on the length scale of d which grows in width as δz ∝ ω 2/3 [14,63]. This broadening is seen in figure 2 as the gap between the turning point located at z c and the horizon at z = 0, and also as the region between the solid black and dashed green curves in figure 5.
We identify five distinct points a) through e) as indicated in figure 6, and will apply our method at these points for each choice of branch cut, starting with the +iπ cut, although our method can just as well be applied for the −iπ choice (as will later be shown) and to the ω < 0 half plane. Each steepest descent contour diagram is created in the s integration plane for a particular set of control parameters (z, ω). To ensure convergence of the integral the real s contour must be deformed into paths of steepest ascent/descent starting and ending in asymptotic valleys at infinity V 1 and V 3 respectively where Re[λf (s, z, ω)] → +∞, |s| → +∞, passing through saddles s j , j = 1, 2, 3 of the phase and satisfying The steepest descent contours may take excursions to and from intermediate valleys, for example, V 2 or its copies on different Riemann sheets. The arguments of s that determine the asymptotic valleys will depend on which choice of cut is taken for the logarithm, but copies on different sheets will always have values of arg(s) separated by 2π. Different subsets of saddles s j can contribute at different values of (z, ω). This is as a result of the topology of the steepest descent paths changing at a Stokes's line in the (z, ω) plane. The steepest paths emerging from two (or more) saddles i = j connect in the complex s-plane when If the number of saddlepoints contributing to the asymptotic expansion of the integral changes as a Stokes line is crossed this is termed a Stokes's phenomenon, see for example [5].

+iπ branch cut
We shall first focus on the region characterized by points e) in figure 6. The steepest descent contour diagram of Eq. (13) at these points is given by the left hand plot in figure 7. The presence of the cut complicates progress in the s-plane. The dashed black line indicates the +iπ branch cut, and at first glance there seems to be no obvious choice of contour (solid black line) starting in V 1 , passing through any of the saddles (black points), and ending in V 3 .
CPF [63] proceeded by effectively removing the logarithmic term from the phase (which requires ω to be small) for the purposes of saddlepoint analysis such that the analysis reduces to that of the Airy function modified by the logarithmic branch cut. They allowed their deformed steepest descent contours to snag on the branch cut and expanded asymptotically around that loop contour. This loop contour is not along a path of steepest descent, but is approximately evaluated in terms of a complex gamma function by the "dominated convergence theorem" (k c → ∞), which to zeroth order ignores nonlinear dispersive effects.
Here, we instead continue to follow the steepest paths, even as they encounter the branch cut and flow onto adjacent Riemann sheets. The result is a calculation that then relies just on simple expansions around saddlepoints including, where required, on the non-principal sheet. This has a potential for an easier physical interpretation than the loop contour around the cut, and additionally does not neglect any near horizon behavior since the nonlinearity in the dispersion is taken into account. We discuss the relative differences of the approach of CPF [63] and this paper further in the discussion below.
To that end, we follow Stone et al [43], and apply the exponential transformation to the integrand of Eq. (13) to give This unfolds the infinite number of Riemann sheets arising from the logarithm in the complex s-space, whose boundaries are defined by the choice of +iπ branch cut, into a single complex w-space. It also removes the need to consider the pole contribution in the original integral for ω = 0 (when ω = 0, this transformation generates a third saddlepoint at ω → −∞ equivalent to the pole). Each sheet is mapped to a horizontally stacked semi-infinite strip of height 2πi in the w-plane, starting with this choice of branch cut at w = +iπ. As a consequence, the w-plane contains an infinite number of periodically stacked valleys V j,n and saddlepoints w j,n = w j,0 + 2nπi, (21) where j = 1, 2, 3, and n ∈ Z denotes the Riemann sheet on which the valleys or saddles sit (n = 0 being the principal sheet for the +iπ choice of branch cut). There is a basis of 3 saddlepoints on the principal sheet, with copies equally spaced out at a complex distance 2πi on each mapped sheet. The w-plane representation is shown in the right hand plot in figure 7. The solid black lines in the w-plane denote the images of the steepest descent contours. The horizontal black dashed lines denote the mappings of the (now no longer) +iπ branch cut. The blue regions denote the asymptotic valleys of convergence, and the red regions are the asymptotic hills where the integration along a contour would diverge. In w-space the periodically repeating valleys V j,n all lie along the positive Re(w) side as w → ∞. The arrows in both plots of figure 7 indicate the direction of travel along the now continuously deformed steepest path staring from V 1,0 in the principal sheet, passing through a subset of the saddles w j,n , on different Riemann sheets if needs be, before ending back at V 3,0 back on the principal sheet.
In w-space starting at V 1,0 within the principal sheet, the contour intersects the first saddle w 1,0 and runs off into the V 2,0 valley "above it", along the bottom of the branch cut at w = +iπ. The deformed contour then re-emerges from V 2,0 intersects a second saddle w 2,0 , before leaving the principal sheet to encounter w 2,−1 which is a copy of the saddle w 2,0 at a point 2πi vertically below on the next lowest sheet. The contour then turns by a right angle (indicating Stokes's phenomenon) before passing into V 2,−1 which is a copy of V 2,0 on the lower side of border with the principal sheet. The last component of the contour re-emerges from V 2,−1 and passes through saddle w 3,0 before finally running into the asymptotic valley V 3,0 on the principal sheet. The contours are considerably easier to follow in w-space than s-space, and it is seen that all asymptotic contributions to the integral arise from saddlepoints (or the pole outside the exponential), rather than loops around branch cuts, allowing for a (local) application of catastrophe theory (the topological theory underlying the coalesence of stationary points) to understand the properties of the integral.  figure 6, together with the choice of +iπ branch cut (dashed black line). The black dots denote the saddlepoints s j of the phase, given by Eq. (16), and the arrows indicate the contributing (converging) steepest descent contours. Right: Equivalent steepest descent contour diagram in the transformed w-space. The horizontal dashed lines denote the mappings of the (now no longer) +iπ branch cut, and the black dots are the unfolded saddlepoints w j,n . In both the left and right panels the larger arrows indicate contributing contours on the principal sheet, and the smaller arrows indicate those contributing on adjacent sheets. It is clear that it is easier to follow the contributing contour in w-space than in s-space.

Asymptotic contributions
In terms of asymptotic contributions from saddlepoints in the original s and transformed w-plane, we have the correspondence: where s 2,−1 is the image of saddle s 2 on the next lowest Riemann sheet in the s-plane. The asymptotic contribution from the expansion about each saddlepoint w j,n along the doubly infinite steepest decent contour that passes through it takes the form [8] ψ (j,n) (z, ω) = where, with f j,n = f (w j,n , z, ω), N j,n is the number of terms taken in the asymptotic series expansion for the corresponding w j,n saddlepoint, Γ denotes the gamma function, g(w) = 1/(2π) which does not actually depend on w and is constant in our case, and q = 0 or 1, depending on the direction of traversal of the contour relative to the computed T (j,n) r . When the expansion is undertaken over a semi-infinite contour, as occurs at a Stokes's phenomenon, additional terms at half powers of the large asymptotic parameter are required, see [78,79].
The terms in the expansion may be computed via the residue integral representation of Eq. (24) or via the Lagrange reversion method [80]. The first couple of terms in a doubly infinite contour expansion over w j,n are: The expressions for f j,n are algebraically complicated and will not be written out here but they, and hence the associated saddlepoint expansions, are valid for a general range of ω. However, for the purposes of comparison with [63], we only need their analytical form in the small ω limit (this is also the regime where the coupling between the positive and negative wavenumber solutions that gives the Hawking effect is strongest). In order to demonstrate our method we focus our attention on points e) which lie below the caustic in figure 6, although it can and will be later applied to the remaining points (it could also be applied for ω < 0). Points e) correspond to saddlepoint diagrams equivalent to that of figure 7. At such points, the principal sheet contributes three real (in s-space) saddles, describing three wave interference for ω > 0 inside of the horizon and below the caustic, like in the Pearcey function. However, here the outer two saddles in the left-hand plot of figure 7 give contributions that dominate the third, leading to the Airy-like interference pattern observed in the ω > 0 region of the ψ +iπ plots in figure  5. Furthermore, this approach also reveals the presence of a contribution from a fourth saddle, w 2,−1 . In this way we find that for the +iπ branch cut at points e), the formal asymptotic expansion ψ (e) +iπ (z, ω) takes the form The contribution from w 2,−1 is exponentially smaller than that from w 2,0 by a factor e −2λπω (ω > 0). Such real exponential factors are associated with pair production [46,54,63]. The factor of 1/2 is due to the presence of the Stokes's phenomenon which may be inferred from the contour intersecting this saddle making a sudden sharp "dogleg" turn as it encounters the saddle (see the right hand panel in figure 7), a characteristic signature of a Stokes's phenomenon.
Due to the fact that the T (j,n) r (z, ω) effectively only depend on derivatives of f (w, z, ω) at w j,n , it is easy to see that for each r, Taking into account the relative sense of the traversal of the contours over w 2,0 and w 2,−1 we find that Eq. (27) simplifies at leading order to (29) From Eqns. (10) and (14) we see that in the notation of [63], their ω/κ factors are equivalent to our λω factors. In order to compare against the results of CPF we observe that we need to consider the small ω regime: taking the limit ω → 0 + of f j,n and T (j,n) 0 , we find that the factors that contribute to Eq. (29) can be written as Up to overall prefactors (see also the next paragraph), Eqns. (30) and (32) are consistent with the results given in equations 62 and 63 of [63] respectively, with the identification of CPF's ∆(z) = |z| 3/2 . The removal in [63] of the log term from the exponent before undertaking a steepest descent approach assumes both |z| 1 and |z| ω 2/3 [46,63], and is consistent with the small ω approximation made here. However, our retention of the logarithm in the exponent increases the range of validity of the results in ω (albeit at the expense of algebraic complexity). A consequence of removing the logarithm from the saddlepoint exponent in [63] is that although this gives an accurate asymptotic approximation for large enough z, it does not capture the presence of a separate horizon and cusp caustic.
As we now explain, our Eq. (31), multiplied by the (1 + e −2λπω /2) factor from Eq. (29), is equivalent to CPF's more complicated looking third contribution which is given by equation 64 in [63]. In that paper equation 64 does not come from a saddlepoint contribution but instead from a loop contour around the cut along -iR in the s-plane (the left image in their figure 4). The authors then apply the dominated convergence theorem which requires that k c → ∞ (equivalent to our λ → ∞) and so to lowest order ignores the cubic term in the integral exponent. After notational translations, equation 64 in [63] contains terms involving sinh(λωπ)Γ(−iλω). Their result expands on the Stokes line of the gamma function according to equation 3.4 of [81] for z = −iλω, with λω > 0 as: Hence for large λ the leading order result of equation 64 in [63] yields Eq. (31). The additional subdominant contribution in the second line of Eq. (33), +(1/2)e −2πiz = +(1/2)e −2λπω , corresponds to the contribution of the fourth saddle w 2,−1 which, since f (2,0) − f (2,−1) = 2πi, combines to give the prefactor (1 + e −2λπω /2) in the overall w 2,0 term in Eq. (29). From this, it can be seen that the pure-saddlepoint approach not only incorporates the cubic terms in the exponent but also avoids the need for complex gamma functions (and the complexity of the correct analytical representation of them on their Stokes lines) whether for small or more general values of ω.
Physically speaking, and for values of z below the cusp, the contributions w 1,0 [Eq. (30)] and w 3,0 [Eq. (32)] correspond to the two waves φ in −k and φ in k respectively, in the right plot of figure 2 (these waves are responsible for the Airy-like interference in the wave plots shown in figure 5). The combination of w 2,0 and w 2,−1 [Eq. (31) multiplied by the (1 + 1/2e −2λπω ) factor from Eq. (29)] generates the φ out −k wave in the same region (responsible for the step function which is modulated by the Airy-like interference). The direction of travel of these three waves can be confirmed by realizing that the spatially dependent parts of the phases in Eqns. (30)-(32) are given by the WKB result z 0 k dz , so that k(z) can be obtained by differentiation, and the group velocity found from v gp = (dk/dω) −1 . For the waves in Eqns. (30) and (32) this gives v gp = 2κd|z|/λ confirming that they are right moving. For the wave in Eq. (31) it is v gp = −κd|z|/λ showing that it is left moving. Interestingly, the group velocity does not depend on the 'Airy factors' ±(2/3)iλ|z| 3/2 as these do not depend on ω. Figure 8: Summary of the saddlepoint contributions at points a) through e) from figure  6, for the +iπ branch cut. The upper plots are the complex s-space contour diagrams, while the lower plots are the corresponding w-space ones. We use the notation s i+j (i = j) to denote when multiple saddlepoints s i and s j coalesce (which happens at a caustic) and become equal to one another. We also adopt this notation for the w saddles.

+iπ branch cut contributions
We can proceed in the same way for each point in the (z, ω > 0) half-plane. The qualitative steepest paths and the associated contributing saddles at these locations are displayed in figure 8. Along the caustic at point d), as expected, two of the saddles from e) (figure 7 or the rightmost column in figure 8) have coalesced into one, w 1+2,0 , and contribute whilst simultaneously undergoing a Stokes's phenomenon with an exponentially subdominant pair of coalesced copies w 1+2,−1 , together with a simple real saddle w 3,0 with a purely imaginary phase.
For values of z that lie between the caustic and the horizon, corresponding to points c), there is one real saddle w 3,0 and one complex saddle w 1,0 undergoing a continuous Stokes's phenomenon with the contributing subdominant saddle w 2,−1 . The latter saddle lies outside of the principle Riemann sheet, and all together the contributing saddles yield ψ (c) +iπ . By a continuous Stokes's phenomenon we mean it occurs for a range of z as opposed to at a single value of z.
At the horizon z = 0, which is point b), w 1,0 now undergoes an instantaneous Stokes's phenomenon (i.e. only at the single point z = 0) inside the principle sheet with the subdominant contributing saddle w 3,0 . Finally, at points a) we see the real saddle no longer contributes and only a single complex saddle w 1,0 remains (corresponding to φ out ↓ in the right plot of figure 2). Consideration of the steepest contours in figure 8 shows that there is a change in the number of contributions across the event horizon at z = 0 (w 1,0 contributes on both sides, but w 3,0 only contributes for z < 0). Hence the event horizon is a Stokes line for real ω.
This is consistent with WKB analysis in a single complex dimension, where Stokes lines sprout from turning points, of which the caustic here is a higher dimensional version. The connection between certain types of black holes, horizons, pair-production, and Stokes's phenomenon has also been studied in various other contexts [82,83,84,85]. In our case the Stokes surface differs: it intersects real 2-parameter space in a line corresponding with the horizon.

−iπ branch cut contributions
We can proceed in the same way for the −iπ choice of branch cut. For both s-and w-space, the contours and saddles in figure 9 are exactly the same as those in figure 8. However, the different location of branch cut forces the starting valley V 1,0 to be shifted to its copy V 1,−1 , a distance 2πi below in the w integration plane. In other words, the principle sheet for the −iπ cut has been shifted and differs from that of the +iπ cut, as can be seen by comparing figures 8 and 9. At points e) the two saddles w 1,−1 and w 3,0 contribute (again corresponding to φ in −k and φ in k , respectively), in contrast to the four contributions from the analogous (z, ω) point for the +iπ branch cut. This describes the wave interference we see in the right plot of figure 5. At the caustic d), two of these saddles w 1,−1 and w 2,−1 coalesce, so that there is one double saddle and one simple saddle contribution. Above the caustic at c) a Stokes's phenomenon is continuously occurring between w 1,−1 and w 2,−1 , (w 1,−1 is subdominant to w 2,−1 ). This persists until the horizon at z = 0 is the reached, at which point a double Stokes's phenomenon takes place at point b) (two dog leg turns on the steepest path from w 1,−1 to w 2,−1 to w 3,0 ), with w 3,0 dominant, w 2,−1 subdominant and w 1,−1 sub-subdominant. Finally, beyond the horizon at points a), the contour over w 1,−1 (φ out ↓ ) passes onto the next lowest sheet, running into and out of V 3,−2 before encountering w 3,−1 (φ out k ) turns through a right angle before running back up to the principal sheet, passing over w 3,0 before running into V 3,0 . From this we observe that w 3,−1 is always undergoing a Stokes's phenomenon above the horizon for z > 0, and is subdominant when compared to w 3,0 .
We make the following remarks: First, the two choices of branches +iπ and −iπ display complementary contributions from sub-subdominant saddles located outside their respective principal Riemann sheets. Whenever the +iπ contour diagrams have contributing saddles outside the principal sheet the −iπ diagrams do not, and vice versa. Second, at points a), aside from the a single real wave, we have two additional waves: the evanescent wave from within the principal sheet and the sub-subdominant saddle from an adjacent one. This differs from the results of CPF [63], where they only find there to be a single real saddle. This is because of their exclusion of the cubic term from the phase (k c → ∞), prior to applying the dominated convergence theorem (see equation 66 in [63]). The two approaches would agree if exponentially subdominant asymptotic contributions were to be neglected in the presence of more dominant ones.
Third, along the horizon at z = 0, a double Stokes's phenomenon takes place, as three saddlepoints are simultaneously joined by a single steepest path. This is also an   indication of the potential for a higher order Stokes's phenomenon [86], which would lead to additional interesting behaviour in the non-physical complex (z, ω) space. The horizon is indeed a part of a Stokes set (intersecting the already-identified real (z, ω) Stokes surface), no matter the choice of cut. The overall qualitative behavior of the contributing saddlepoints (waves) is summarized in figure 10 for both choices of cut, and describes the wave behavior observed in figures 5. The number of contributing real and complex saddlepoints together with their relative dominance throughout the (z, ω) plane is shown. Red lines and text denotes contributing saddles which are a part of the Stokes set, whether it be from within or outside of the principle Riemann sheet. The dominance of the saddles is denoted by the " > " signs. These plots can be compared against figure 4 for the Pearcey function. Clearly, the event horizon catastrophe has considerably more structure.

Numerical Comparison and Validity of Asymptotics
In figures 11 and 12 we demonstrate that our exponential approach yields an accurate approximation by plotting the asymptotic expressions against the exact wave functions ψ ±iπ obtained by numerically solving the integral in Eq. (13) for the two choices of +iπ which is only valid at the horizon, and we do not label ψ (d) +iπ since point d) corresponds to the caustic where the asymptotic approximation diverges. As expected, the asymptotic expansions provide an (exponentially) good approximation to the exact result except near the caustic at z c = −3/(10 2/3 ). The slight disagreement in the asymptotics at the horizon will vanish for large λ → ∞. The saddles w j,n that the steepest paths encounter in each region are also indicated, together with their relative dominance (denoted by ">"). The saddles w 1,0 and w 3,0 have equal real parts for z below the caustic (denoted by "∼").  (13) but now for the −iπ choice of branch cut. The notation and selected parameters are the same as those used in figure 11, and similarly this represents a constant ω > 0 slice of the righthand plots in figure 5. The zeroth r = 0 order asymptotic approximations a) through e) (solid red lines) and optimally truncated asymptotic approximations (blue dots, see Appendix A) are again plotted for comparison. As expected, the asymptotic expansions provide an (exponentially) good approximation to the exact result except near the caustic at z c = −3/(10 2/3 ) and exactly at the horizon z = 0. The saddles w 1,−1 and w 3,0 have equal real parts for z below the caustic (denoted by "∼").
terms r = 0 in the asymptotic expansions [see Eq. (23)] and yet find an excellent match to the exact results (sufficiently far from the caustic). This is despite the fact that our 'large' parameter λ is only of order unity. The blow up close to the caustic could be fixed by using a uniform approximation [68,87]. It is also noteworthy that the asymptotic approximations for the different spatial regions match together so well at their respective borders. The only slight discrepancy is at the horizon between the asymptotic values of ψ +iπ . We find that when we use larger values of λ this difference vanishes as expected so that ψ Eq. (23) is a diverging series, so there is an optimal series truncation in r which can be made for each saddlepoint contribution, giving the approximation of lowest possible error. Taking terms higher than this value of r will actually decrease the validity of the approximation and increase the error. Although we find that the lowest order r = 0 truncation already gives a good visual match, we have also evaluated the higher terms and these are included as the blue dots in figures 11 and 12. The details of the calculations of the higher order terms are given in Appendix A of this work.

Concluding remarks
In this paper we take a 'catastrophe theory' approach to horizons motivated by the observation that a nonlinear dispersion relation causes the solutions of Hamilton's equations to undergo a broken pitchfork bifurcation near the horizon. Pitchfork bifurcations are specified by two control parameters which in our case are the lab frame frequency ω and the position coordinate z.
According to catastrophe theory, the universal structurally stable relationship between these control parameters gives a cusp shape z ∝ −ω 2/3 which defines the location of a caustic z c where waves coalesce. However, whereas cusp caustics are usually dressed by the Pearcey function wave pattern, the event horizon bifurcation gives rise to a novel form of wave pattern described by an Airy function modified by a logarithmic term we call the log-Airy function.
Some familiar properties remain such as self-similar scaling and we use this to identify a classical limit with a linear dispersion. Furthermore, like the Pearcey function there is a Stokes set that occurs outside the cusp, although in the log-Airy case it is flattened into a straight line in the (ω, z) plane which coincides with the event horizon. Except for the special case of ω = 0, the caustic at z c [point d) in figure 6] and event horizon at z = 0 [point b) in figure 6] do not sit at the same location: the caustic lies downstream behind the horizon and the shape of the caustic implies that the spatial gap between them grows as ω 2/3 . This scaling has been pointed out before on the basis of the behaviour of the classical solutions [63], and the connection between caustics and horizons was previously studied in a different way in the context of water waves [16,49]. However, the knowledge that it is a universal prediction of catastrophe theory and corresponds to the zone between a Stokes line and caustic adds to our understanding of the notion of a broadened horizon. On the other hand, this challenges us to generalize wave catastrophe theory to include logarithmic terms that ultimately arise from particle creation in quantum field theory [47,88,89].
The log-Airy function has previously been analyzed by CPF in [63]. However, our treatment differs in some important respects. In CPF [63] the expansion in the region characterized by point e) (see figure 6) is based on two saddlepoints plus a loop contour around a branch cut. The effective removal of the logarithm by CPF from the phase significantly simplifies expressions for their two saddlepoint contributions, but restricts the validity of their results to |z| d and |z| d (ω/κ) 2/3 (units restored), see equations 34 and 57 in [63] and equations 11a and 11b in [46]. Their two saddlepoint contributions are equivalent to our Eqns. (30) and (32), which were obtained by performing a small ω expansion to our large λ asymptotic result, Eq. (29). In fact keeping the logarithm in the exponent is necessary to understand the role the latter plays as a Stokes line across which the number of contributing saddlepoints (and so waves) changes. The result of their loop contribution (corresponding to our Eq. (31) in the large λ limit) is a complicated cut expansion involving the gamma function on its complex Stokes line. Their expression approximates k c → ∞, ignoring nonlinear effects, and obscures the underlying simplicity of the contributing subdominant copies of saddlepoints on adjacent Riemann sheets.
Although we only explicitly focused on an analytic description for points e) [Eq. (29)] for the +iπ choice of branch cut, our method was applied to the remaining points for both choices of cut and for a small constant value of ω > 0 in order to show the validity of our approach. This is shown in figures 11 and 12 for both the zeroth order analytic approximations [obtained via Eq. (23)] and for the optimally truncated numerical asymptotics (described in Appendix A). Again, this approach could just as well be applied for larger ω and for ω < 0.
Our asymptotic contributions still diverge at the caustic as expected, but a uniform approximation [68,87] could be employed locally to deal with this. Due to the divergence being attributable to a coalescence of two saddlepoints, this would take the form of an Airy function and its derivative. The caustic at which this occurs in the (z, ω) plane is one dimensional and hence even though there is a logarithm perturbing the polynomial in the exponent of the integral, the local behaviour across this cusp is still structurally stable and so falls within the realm of catastrophe theory.
Our large parameter λ given in Eq. (15) is a constant that does not depend on position, unlike ∆(z) defined in equation 34 of [63] which vanishes at the horizon. Combined with our transformation to exponential coordinates, this allows the near horizon behaviour to be examined with exponential accuracy, including the elucidation of new sub-subdominant contributions in regions below the caustic. Furthermore, the present catastrophe motivated approach could allow tight bounds to be put on corrections to the Hawking spectrum of emitted particles due to nonlinear dispersive effects, particularly if ω is not small. In particular, CPF [63] show how to combine the wave functions ψ ±iπ to obtain the Bogoliubov coefficients that directly give the Hawking production rate. This will be pursued in future work. a welcoming and nourishing environment in Bristol. We fondly remember our first encounters with wave singularities not just on paper but in real life in the form of expeditions to see the Severn bore. Chasing across the Gloucestershire countryside in Michael's car, sometimes at night, to an ideal viewing spot we recall the low roar of the wave as it approaches, which is especially intimidating after dark, accompanied by reversal of the river's flow direction, change in depth and sometimes wet feet. It provides observers with an impressive and tangible demonstration of Nature's power, not to mention the universality of the Airy function.