## Abstract

This paper presents a method of analysing wave-field dispersion relations in which Bessel functions of imaginary order occur. Such dispersion relations arise in applied studies in oceanography and astronomy, for example. The method involves the asymptotic theory developed by Dunster in 1990, and leads to simple analytical approximations containing only trigonometric and exponential functions. Comparisons with accurate numerical calculations show that the resulting approximations to the dispersion relation are highly accurate. In particular, the approximations are powerful enough to reveal the fine structure in the dispersion relation and so identify different wave regimes corresponding to different balances of physical processes. Details of the method are presented for the fluid-dynamical problem that stimulated this analysis, namely the dynamics of an internal ocean wave in the presence of an aerated surface layer; the method identifies and gives different approximations for the subcritical, supercritical and critical regimes. The method is potentially useful in a wide range of problems in wave theory and stability theory. A mathematical theme of the paper is that of the removable singularity.

## 1. Introduction

The aim of this study was to present a hierarchy of analytical approximations to a type of dispersion relation that can occur in applied studies of wave motion and that is complicated by the presence of Bessel functions of imaginary order. A guiding principle of the paper is that the final approximations should involve only elementary functions, and be easy to analyse by research-workers who would not themselves wish to explore the intricacies of the special functions of mathematical physics. One such intricacy is that both the argument and order of the Bessel functions may be non-analytic functions of frequency and wavenumber, typically involving square roots, even though the dispersion relation taken as a whole is analytic. Thus the exact dispersion relation is a hard-to-interpret combination of Bessel functions of real or complex orders, in which unexpected combinations and cancellations of terms occur in some regimes but not in others. Nevertheless, it will be shown that a hierarchy of approximations can be derived that involve only elementary functions and that yet have extremely high accuracy. The asymptotic formulae of Dunster (1990) for Bessel functions of imaginary order are invaluable for this purpose. A hierarchy of approximations to the dispersion relation is worth having, not only because it makes possible a simple description of parametric dependence, but also because it identifies regimes in which different physical effects are dominant or negligible. The high accuracy of the approximations makes them well-suited to computational work, but the main thrust of the paper is that very simple and accurate analytic approximations are available where they might not be expected.

Dispersion relations of the type considered here arise particularly in studies of wave propagation in stratified media—for example, in oceanography (Brekhovskikh & Godin 1990) and in astronomy (Diaz *et al.* 2007). Linear or exponential profiles often lead to differential equations with solutions expressible in terms of confluent hypergeometric functions, within which Bessel functions form a large sub-class.

The paper is organized as follows. Section 2 summarizes the required theory for a non-trivial example of wave propagation, that of an internal ocean wave in the presence of an aerated surface layer; this physical problem, as presented in Grimshaw *et al.* (2008), was the stimulus for the present work. A critical frequency is identified, defined mathematically by the non-analytic transition of the order of a Bessel function from imaginary to real values and defined physically by the transition of the wave structure from that of a modified internal wave to that of a modified aerated surface wave. This frequency separates the subcritical regime analysed in §3 from the supercritical regime analysed in §4. Close to the critical frequency, i.e. in the critical regime, the coupling of the two types of wave is at its strongest and gives maximally interacting waves; the dispersion relation for this regime is analysed in §5. Mathematically, the critical regime corresponds in the (wavenumber, frequency) plane to the neighbourhood of a removable singularity in the dispersion relation. Two principal results in the paper are the ‘main subcritical approximation’, also called the Dunster–Taylor approximation, derived in §3, and the ‘main supercritical approximation’, which is a Debye approximation, derived in §4. Further approximations are possible, some with a narrow range of validity; these are indicated at various points in the paper. Conclusions are presented in §6.

## 2. Internal waves with an aerated surface layer

This section gives the governing equations for internal waves coupled to an aerated surface layer in which bubbles are present, and derives the dispersion relation for such waves (Grimshaw *et al.* 2008). An intrinsic frequency that distinguishes different regimes and makes possible a physically lucid set of scalings and non-dimensional variables is identified.

### (a) Internal waves

We begin by specifying the characteristics of a stratified layer of incompressible liquid without aeration. The layer is −*H*≤*z*≤0, with rigid lower surface *z*=−*H* and free upper surface *z*=0. In the motionless base state, the density of the liquid is *ρ*_{l0}(*z*), in which the subscript l denotes non-aerated liquid and the subscript 0 denotes the base state. Thus the buoyancy frequency *N*_{l} is given by *N*^{2}_{l}=−*gρ*_{l0z}/*ρ*_{l0}, with *g* denoting the acceleration due to gravity and the subscript *z* denoting differentiation with respect to *z*. We assume that *N*_{l} is constant, possibly zero, so that
2.1where *ρ*_{a} is the density at the upper surface *z*=0. The pressure *p*_{l0}(*z*) of the liquid in the base state satisfies the hydrostatic relation *p*_{l0z}=−*gρ*_{l0}, so that
2.2where *p*_{a} is atmospheric pressure.

Internal waves perturb the motionless base state of non-aerated liquid to give a vertical component of velocity , which we take in the modal form
2.3where the real part is understood. Thus we consider two-dimensional flow in which *x* is the horizontal coordinate, *k* is the horizontal wavenumber, *t* is time and *ω* is the frequency. Standard Boussinesq theory (Lighthill 1978; Gill 1982) gives the Sturm–Liouville problem
2.4Here the boundary condition *w*(0)=0 is the rigid-lid approximation for a free surface. For given values of *N*_{l} and *H*, the problem defined by equations (2.4) determines the dispersion relation *D*(*ω*,*k*)=0 characterizing the modes.

### (b) Aerated surface-layer waves

We now consider the effect of an aerated surface layer. In the motionless base state, the degree of aeration as a function of depth is specified by the void fraction *α*_{g0}(*z*), defined so that the density of the aerated liquid is *ρ*_{l0}(*z*)(1−*α*_{g0}(*z*)). The subscript g denotes gas, and it is assumed that *α*_{g0}≪1. For typical oceanic conditions, the main effect of aeration on the buoyancy frequency is that the quantity *ρ*_{l0z}/*ρ*_{l0}, representing the reciprocal of a vertical length-scale of stratification in the pure liquid, must be replaced by its aerated counterpart *ρ*_{l0z}/*ρ*_{l0}−*α*_{g0z}. Hence internal waves in the presence of an aerated surface layer are typically still governed by equations (2.4), provided that the first equation is replaced by
2.5where *N*_{eff}(*z*) is the effective buoyancy frequency defined by
2.6We shall assume the void fraction varies exponentially with depth according to , where *α*_{a} is the void fraction at the upper surface *z*=0, and *h*_{g} is the skin-depth of the aerated layer. Thus
2.7

### (c) Scaling and non-dimensional variables

As indicated in §2*a*, we take the buoyancy frequency *N*_{l} to be a constant. The aerated wave problem posed in §2*b* contains the length-scale *h*_{g}, and the form of equation (2.7) shows that the problem has an intrinsic frequency scale *ω*_{g} defined by
2.8Thus we may introduce the non-dimensional variables
2.9The coupling in the problem is represented by the frequency ratio *N*_{l}/*ω*_{g}. This determines a dimensionless buoyancy frequency *Ω*_{N} defined by
2.10Other dimensionless constants we shall use are the small quantities *δ* and *ϵ* defined by
2.11Thus *δ* is the fraction of the total fluid depth occupied by the aerated surface layer. Although, at first sight, *ϵ* is too small to be a useful quantity, it is convenient in intermediate formulae.

The dimensionless buoyancy frequency *Ω*_{N} may take a wide range of values. A typical oceanographic value is about 0.2, that is *N*_{l}≃0.2*ω*_{g}. This could occur with, for example, *h*_{g}=4 m and *α*_{a}=10^{−3}, for which *ω*_{g}=0.07 s^{−1}; a typical value *N*_{l}=0.015 s^{−1} then gives *Ω*_{N}=0.21. The dimensionless frequency *Ω* is such that the range *Ω*<*Ω*_{N} corresponds to *ω*<*N*_{l}, and the range *Ω*_{N}<*Ω*<1 corresponds to *N*_{l}<*ω*<*ω*_{g}. We shall call the former range the subcritical regime, and the latter range the supercritical regime. The critical frequency *Ω*=*Ω*_{N} is simply the buoyancy frequency *ω*=*N*_{l} in the absence of aeration. The limiting case *Ω*_{N}=0 arising from *N*_{l}=0 is also of interest, as representing internal waves in the aerated layer when there is no background stratification; these are aerated modes.

It is convenient to express depths in terms of the variable defined by
2.12This converts the aerated wave problem based on equations (2.5) and (2.7) to the Bessel-equation form
2.13Here and
2.14The general solution of (2.13)_{1} is an arbitrary linear combination of the Bessel functions , or equivalently, when *ν* is not an integer, of . Application of the boundary conditions in (2.13) gives the dispersion relation in the equivalent forms
2.15or
2.16This is the dispersion relation in Grimshaw *et al.* (2008), expressed in variables that are most suited for relating the various asymptotic regimes to the underlying physical processes. In the subcritical regime *Ω*<*Ω*_{N}, that is *ω*<*N*_{l}, the quantity *ν*^{2} defined by equation (2.14) is negative, and so *ν* in equations (2.15) and (2.16) is imaginary. Thus, Bessel functions of imaginary order are unavoidable. Simple but highly accurate approximations to these functions are given in Dunster (1990), from which a selection is presented in Olver *et al.* (2009). For the basic theory, see Olver (1974).

Series expansion of *J*_{±ν}, followed by cancellation of (*K*/*Ω*)^{±ν} and division by *ν*, casts equation (2.16) in a form containing only even powers of *ν*. Thus, the apparent singularities are removable; hence despite initial appearances the dispersion relation is analytic in *Ω* and *K*, and can be written in a series form containing neither square roots nor *i*. Although this series form shows that the dispersion relation is analytic, it is not easily related to standard formulae, and at large argument is numerically ill-conditioned. For our purposes, it is better to use equations (2.15) and (2.16) as they stand, taking as a starting point whichever is more convenient in a particular calculation.

### (d) Scope of the present work

The opaqueness of the dispersion relation (2.15) or (2.16) is in striking contrast to the simplicity of the dispersion relation for internal waves in a uniformly stratified layer without aeration. Yet on physical grounds, in the subcritical regime *Ω*<*Ω*_{N}, the dispersion relation (2.15) or (2.16) must be a modification of this simpler form, with the degree of modification increasing as the critical frequency *Ω*_{N} is approached from below. Equally, in the supercritical regime *Ω*_{N}<*Ω*<1, the dispersion relation must be a modification of the dispersion relation for an aerated layer in the absence of background stratification, with the degree of modification increasing as the critical frequency is approached from above. These remarks follow from the basic property of internal waves that they have an upper cut-off frequency; this frequency is *ω*=*N*_{l}, i.e. *Ω*=*Ω*_{N}, in the absence of aeration, and *ω*=*ω*_{g}, i.e. *Ω*=1, in the absence of background stratification (Lighthill 1978; Gill 1982).

The present work derives simple analytical approximations to the dispersion relation (2.15) or (2.16) for the different regimes. These approximations have physical interpretations based on the type of reasoning just given, and make explicit the way in which the shape of the dispersion curves in the (*Ω*,*K*) plane depends parametrically on *ϵ* and *Ω*_{N}. Comparisons with accurate numerical computations from equation (2.15) or (2.16) are presented to show that the approximations are extraordinarily accurate. The numerical computations, as opposed to the analytical approximations, require software for numerical evaluation of Bessel functions of complex order. In the numerical computations, a practical point is the extreme smallness of *ϵ*; for example, if , then *ϵ*=e^{−25}≃1.4×10^{−11}. This smallness is not present in the final results, which depend not on *ϵ* directly but rather on , i.e. *δ*.

## 3. The subcritical regime: modified internal waves

### (a) The exact subcritical dispersion relation

In the subcritical regime *Ω*<*Ω*_{N}, for which the order *ν* of the Bessel functions is imaginary, we take *ν*=i*ν*_{i}, where *ν*_{i} is real and positive. As this paper is concerned only with propagating modes, rather than with evanescent modes, all wavenumbers and frequencies are real, so that the Bessel functions have real argument. We assume this henceforth. With the notation r for real part, and i for imaginary part, we then have
3.1where *J*_{iνi,i} and *Y* _{iνi,i} are real. Two identities are
3.2These formulae enable the dispersion relation (2.15) or (2.16) to be written, after an exact cancellation of the imaginary part, in either of the real forms
3.3and
3.4The basic functions in Dunster (1990) are not *J*_{iνi,r} and *Y* _{iνi,r}, but rather the Dunster functions and defined for real *ν*_{i} by the relations (Olver *et al.* 2009, p. 248)
3.5Identities (3.2) give similar expressions for *Y* _{iνi,r} and *Y* _{iνi,i} in terms of and . In the notation of Dunster (1990), the functions and would be denoted *F*_{iνi} and *G*_{iνi}; we use here the standard later notation of Olver *et al.* (2009). These references have *ν* in place of our *ν*_{i}, and so write the original Bessel functions as *J*_{iν} and *Y* _{iν} rather than as our *J*_{ν} and *Y* _{ν}. For real argument and order, the Dunster functions are real, and we shall use them henceforth. Thus equation (3.3) or (3.4) becomes
3.6This is the basic equation for the subcritical regime *Ω*<*Ω*_{N}, with
3.7

### (b) The main subcritical approximation

#### (i) Dunster functions with argument *K*/*Ω*

Unless *K*/*Ω* is too small, the terms and in equation (3.6) are accurately represented by the Dunster approximations
3.8in which the absence of *ν*_{i} on the right-hand sides is noteworthy. The analysis in Dunster (1990) shows that if *ν*_{i}<1, the approximations (3.8) require *K*/*Ω*>*O*(*ν*_{i}); numerically they are accurate when *K*/*Ω*>2*ν*_{i}. Although the asymptotic theory leading to approximations (3.8) is for ‘large’ *K*/*Ω*, this does not mean that *K*/*Ω* needs to be large compared with 1; much smaller values of *K*/*Ω* are allowed if *ν*_{i} is small. By contrast, if *ν*_{i}>1, the approximations require *K*/*Ω*>*O*(*ν*^{2}_{i}); estimates of the range of validity may be based on *K*/*Ω*>*ν*^{2}_{i}, to within a factor of order 1. This dependence of range of validity on the value of *ν*_{i} needs to be borne in mind when interpreting the unqualified error term applicable to approximations (3.8), which from Olver *et al.* (2009, p. 248) is *O*((*K*/*Ω*)^{−3/2}).

The equations *ν*_{i}=1 and *ν*^{2}_{i}=*K*/*Ω* give two curves
3.9in the (*K*,*Ω*) plane. In the positive half *K*>0 of the strip 0<*Ω*<*Ω*_{N}, the former curve is everywhere to the left of the latter, i.e. *K*_{1}(*Ω*)<*K*_{2}(*Ω*), and so the two curves separate the half-strip into three regions. The first region, 0<*K*<*K*_{1}(*Ω*), corresponds to *ν*_{i}<1, and because, from equation (3.7), *K*/*Ω*>*ν*_{i}/*Ω*_{N} when *Ω*<*Ω*_{N}, it follows that the Dunster approximations (3.8) hold accurately everywhere in this first region if *Ω*_{N} is less than about 0.5. This condition is satisfied in the oceanographic problem we are considering, for which a typical value of *Ω*_{N} is about 0.2. The second region, *K*_{1}(*Ω*)<*K*<*K*_{2}(*Ω*), corresponds to *ν*_{i}>1 and *K*/*Ω*>*ν*^{2}_{i}, so that the Dunster approximations hold in this region too. The third region, *K*>*K*_{2}(*Ω*), corresponds to *ν*_{i}>1 and *K*/*Ω*<*ν*^{2}_{i}, for which the Dunster approximations do not hold.

Combining the above results shows that in the positive half *K*>0 of the subcritical region 0<*Ω*<*Ω*_{N}, the Dunster approximations (3.8) hold where *K*<*K*_{2}(*Ω*), to within a factor of order 1. This gives a very large region in the (*K*,*Ω*) plane, because of the disposition of the coefficients *Ω*_{N} in the definition of *K*_{2}(*Ω*). Moreover, the limiting form for *Ω*≪*Ω*_{N} shows that this region includes a generous wide-angle sector near the origin; for example, if *Ω*_{N}=0.2, this sector is 0<*K*<25*Ω*. The bounding values *K*=*K*_{2}(*Ω*) tend rapidly to infinity as *Ω*_{N} is approached, asymptotically as *K*≃1/{2(*Ω*_{N}−*Ω*)}.

#### (ii) Dunster functions with argument *ϵK*/*Ω*

In the opposite direction from (*i*), unless *K*/*Ω* is large enough to approach 1/*ϵ* the terms and are accurately represented by the Taylor approximations (Olver *et al.* 2009, p. 248)
3.10and
3.11The error terms in these two approximations are both of order (*ϵK*/*Ω*)^{2} if *ν*_{i} is bounded away from zero; but if *ν*_{i}→0 the error term in (3.11) becomes of order because of the logarithmic singularity in *Y* _{0}. For *ν*_{i} up to about 0.6, the function is well represented by the first one or two terms of its Taylor series (Olver *et al.* (2009), §25.8.7 with *z*=−i*ν*_{i})
3.12Here *γ* is Euler's constant 0.5772… and *ζ* is the Riemann zeta function defined by . For general *ν*_{i}, the branch of in is chosen to be continuous with equation (3.12); making this choice is called unwrapping the phase.

#### (iii) Combined use of Dunster and Taylor approximations

Use of approximations (3.8)–(3.12) in equation (3.6), keeping only the first term in equation (3.12) and making use of , gives
3.13This is the main subcritical approximation to the dispersion relation for internal waves modified by an aerated layer. The approximation involves only elementary functions: the original Bessel functions, and all other higher functions, have disappeared. Although the approximation is still transcendental and implicit in *K* and *Ω*, it is far simpler and more intelligible than the dispersion relation in its original form.

#### (iv) Numerical performance

The numerical performance of the main subcritical approximation (3.13) is superb. For the numerical work in this paper, we take a lower water depth *H*=100 m, aerated depth *h*_{g}=4 m, buoyancy frequency *N*_{l}=0.015 s^{−1}, upper surface void fraction *α*_{a}=10^{−3} and gravitational acceleration *g*=9.8 m s^{−2}. Thus *ω*_{g}=(2*gα*_{g}/*h*_{g})^{1/2}=0.07 s^{−1}, and the dimensionless parameters used in the plots are *δ*=*h*_{g}/*H*=1/25 and *Ω*_{N}=*N*_{l}/*ω*_{g}=0.21429. In figure 1*a*, the dashed lines are the first 15 branches of equation (3.13) plotted in the (*K*,*Ω*/*Ω*_{N}) plane for 0<*Ω*/*Ω*_{N}<1.2, and the solid lines are the corresponding exact curves computed from equation (2.15). The curves are indistinguishable almost everywhere. In agreement with the analysis in (i), the exact and approximate curves are visibly different only on the branches furthest to the right in the figure, i.e. at large *K*. As the full dispersion curve has infinitely many branches, accumulating on the *K* axis, the branches are not plotted from where they enter a parabolic-shaped region containing this axis. Thus in plots, the termination of branches at their lowest points, where they meet the upper boundary of the parabolic region, is an artefact; all branches, when fully continued, reach the origin (*K*,*Ω*/*Ω*_{N})=(0,0).

An expanded view of the first 10 branches in figure 1*a* is shown in figure 1*b* for the frequency range 0.99<*Ω*/*Ω*_{N}<1.005. Again, the exact and approximate curves are indistinguishable almost everywhere, the only significant difference occurring on arcs close to *Ω*/*Ω*_{N}=1, as discussed below in relation to figure 1*c*. An interesting feature of figure 1*b* is the wavy-stepped structure of the branches somewhat below *Ω*/*Ω*_{N}=1; the stepped structure represents the interaction of the main internal modes with aerated modes in the surface layer. It is striking that the subcritical approximation (3.13) captures this finely detailed structure so accurately. A similar stepped structure occurs in the theory of elastic waves in a layer, where the structure represents the interaction of longitudinal and transverse elastic waves, and the main parts of the steps approach a ‘grid of bounds’ (Mindlin 1960).

The first branch in figure 1*a* is shown in figure 1*c* with a greatly expanded scale of the *K* axis near *Ω*/*Ω*_{N}=1. The figure shows that the limiting error in the subcritical approximation (3.13) as *Ω*/*Ω*_{N}=1 is approached is no greater than about 0.01 in *K*, i.e. 2 per cent, and about 0.001 in *Ω*, i.e. 0.5 per cent. Thus the failure is ‘graceful’, i.e. there is no drastic divergence hidden in approximation (3.13) near *Ω*/*Ω*_{N}=1. This is a consequence of the underlying analytical smoothness of the exact dispersion relation equation (2.15) or (2.16) when the value *Ω*/*Ω*_{N}=1, i.e. *ν*=0, is crossed, as discussed at the end of §2*c*. That is, the apparent singularity at *ν*=0 is removable.

### (c) The internal wave limit

Approximation (3.13) is part of a hierarchy. Subordinate members of the hierarchy are obtained by approximating equation (3.13) further. For example, because *δ*≪1, a simple approximation to equation (3.13) is *ν*_{i}≃*nπδ* for *n*=1,2,…, i.e.
3.14Equivalently, (*Ω*/*Ω*_{N})^{2}=*K*^{2}/{*K*^{2}+(*nπδ*)^{2}}, or in dimensional variables
3.15This is the dispersion relation for internal waves in the absence of aeration. Its numerical performance as an approximation to the exact dispersion relation for coupled internal and aerated waves is good for the first few branches not too far from the origin (*K*,*Ω*)=(0,0), but it loses accuracy as *n* increases, except in regions approaching the origin. This is demonstrated in figure 1*d*, in which the dashed lines are branches of equation (3.14), and the solid lines are the corresponding exact curves. The loss of accuracy is associated only with the approximation of equation (3.13) by equation (3.14), and so casts no doubt on the use of Dunster's approximations (3.8). The loss of accuracy is caused by non-uniformity in *n* and when equation (3.14) is regarded as the first term in a series expansion in *δ*.

### (d) Further approximations in the subcritical hierarchy

An approximation intermediate between equations (3.13) and (3.14) in the hierarchy is obtained by taking *ν*_{i}=*nπδ*(1+*η*), i.e.
3.16When this is substituted into equation (3.13), and the convention is adopted that logarithmic terms are ‘order one’, the leading term in the Taylor series when *δ*≪1 and *η*≪1 gives
3.17Here the form of *η* displays the non-uniformity in *n* and mentioned above, and in consequence, approximation (3.16) provides only a small numerical improvement over approximation (3.14).

Further down the hierarchy even than equation (3.14) is the tangent approximation at the origin, namely *K*/*Ω*≃*nπδ*/*Ω*_{N}. For small *n*, this gives a good approximation to the dispersion curve up to moderate frequencies and wavenumbers, but in the same way as for approximation (3.14), the range of accuracy decreases as *n* increases.

### (e) Assessment of subcritical approximations

Detailed numerical tests have confirmed the privileged place occupied by equation (3.13) in the hierarchy of approximations for subcritical waves. No simplification of equation (3.13) has been found, which can match its almost perfect accuracy over the whole range *Ω*<*Ω*_{N}, as displayed in figure 1*a*–*c*.

## 4. The supercritical regime: modified aerated waves

### (a) Bessel-function forms of the supercritical dispersion relation

In the supercritical regime *Ω*>*Ω*_{N}, the order *ν* of the Bessel functions in the dispersion relation is real, and either of the exact forms of the dispersion relation given in equations (2.15) and (2.16) may be used directly. Because *ϵ* is so small, the terms *Y* _{ν}(*ϵK*/*Ω*) and *J*_{−ν}(*ϵK*/*Ω*) are dominated almost everywhere in the (*K*,*Ω*) plane by their very large factor *ϵ*^{−ν}. Hence the only way the term *J*_{ν}(*K*/*Ω*)*Y* _{ν}(*ϵK*/*Ω*) in equation (2.15) or the term *J*_{ν}(*K*/*Ω*)*J*_{−ν}(*ϵK*/*Ω*) in equation (2.16) can be balanced by the other term in the equation is for
4.1Therefore, this must be the small-*ϵ* limiting form of the supercritical dispersion relation, except for the few special regions in the (*K*,*Ω*)-plane where the above-mentioned dominance does not occur. These are the regions where *K*/*Ω* is large enough to approach 1/*ϵ*, or where *ν*≪1. Because , the latter occurs when the frequency approaches the critical frequency *Ω*=*Ω*_{N}.

The accuracy of approximation (4.1) was checked numerically. In plots comparable to those in figure 1, but for *Ω*>*Ω*_{N}, the branches of equation (4.1) in the (*K*,*Ω*) plane were indistinguishable from the exact branches except in a small region close to *Ω*=*Ω*_{N}. A blow-up of this region on the first branch is shown in figure 1*c*, with a greatly expanded *K* scale. The figure shows that even in this region, the error in *K* is small, no greater than about 0.01, i.e. 2 per cent.

### (b) The main supercritical approximation

In the spirit of our hierarchical approach, we now seek an approximation to equation (4.1) which does not involve Bessel functions. Provided that *K*/*Ω*>*ν* and *K*/*Ω* is not extremely close to *ν*, the function *J*_{ν}(*K*/*Ω*) is well represented by the Debye approximation (Olver *et al.* 2009)
4.2Use of this approximation in equation (4.1) gives
4.3This is the main supercritical approximation, containing only elementary functions. In plots, its accuracy is not visibly different from that of equation (4.1), i.e. the branches of equation (4.3) are indistinguishable from the exact branches except in a tiny region close to *Ω*=*Ω*_{N}. The main difference, indicated in figure 1*c*, is that as *Ω*=*Ω*_{N} is approached, equation (4.3) gives an error in *K* of about 0.02 (i.e. 4%) compared with 0.01 (i.e. 2%) for equation (4.1).

### (c) A further approximation in the supercritical hierarchy

An approximation subordinate to equation (4.3) is obtained by using not the Debye approximation (4.2) but the simpler, though less accurate, Hankel approximation
4.4This may be contrasted with the corresponding Dunster approximation for in equation (3.8)_{1}, in which *ν*_{i} is absent. Use of approximation (4.4) in equation (4.1) gives
4.5Numerical checks show that the accuracy of this approximation degrades rapidly as *Ω* increases, especially for the higher branches. The useful approximation is equation (4.3).

### (d) Another approach to supercritical approximations

Another approach to deriving supercritical approximations is to use the roots *j*_{ν,n} of *J*_{ν}, where *n*=1,2,…, and write the solutions of equation (4.1) in the form *K*/*Ω*≃*j*_{ν,n}. This does not complete the problem, because *ν* is a function of *K* and *Ω*, and *j*_{ν,n} cannot be written as an explicit formula. A variety of approximations is available for *j*_{ν,n}, for example McMahon's approximation (Olver *et al.* 2009). The resulting equations amount to approximations derivable from equation (4.3).

### (e) The limit of no background stratification

It is of interest to examine the limiting case of no background stratification, i.e. *N*_{l}=0, within the framework of this paper. This is the aerated-mode limit. We then have *Ω*_{N}=0, i.e. *ν*=*K* from equation (2.14), and our starting point is equation (2.15) or (2.16) with *ν*=*K* everywhere. There is now no subcritical regime *Ω*<*Ω*_{N}, and the supercritical formulae we have just derived all apply for *Ω*>0 with *ν*=*K*. The approach of §4*d* becomes practical in this case, because it gives *Ω*≃*K*/*j*_{K,n}, in which the right-hand side depends only on *K*. Approximations when *ν*=*K* lose accuracy as the value *Ω*=0 is approached, just as we found above for arbitrary *ν* when the value *Ω*=*Ω*_{N} is approached.

## 5. The critical regime: strong interaction

### (a) Wavenumbers at the critical frequency

In the absence of aeration, the wavenumber of an internal wave tends to infinity as the critical frequency *Ω*=*Ω*_{N} is approached. One effect of aeration is that, on each branch of the dispersion curve, the wavenumber is finite at this frequency. Thus there are infinitely many wavenumbers *K* corresponding to the critical frequency, and we now calculate them.

Since *ν*=0 when *Ω*=*Ω*_{N}, the equation to be solved for *K* is
5.1For small *ϵ*,
5.2In conjunction with the definition *ϵ*=e^{−1/δ}, these expressions give the approximate equation
5.3The omitted term here is of order . Because this is so very small, about 10^{−20} for *δ*=1/25, the error incurred in calculating *K* from equation (5.3) would hardly be detectable in conventional IEEE arithmetic on a computer.

The solutions *K* of equation (5.3) have a simple expansion in powers of *δ*, obtained by starting with the leading-order equation *J*_{0}(*K*/*Ω*_{N})≃0. Hence *K*/*Ω*_{N}≃*j*_{0n} for *n*=1,2,…, where *j*_{01},*j*_{02},… are the zeros of the Bessel function *J*_{0}, namely 2.4048,5.5201,…. The range *n*=1,2,… in *j*_{0n} is assumed henceforth. Higher terms in *δ* are found by taking
5.4and using the Taylor series
5.5and
5.6Substitution into equation (5.3) determines *η*, and hence *K*, as series in *δ*. Identities *J*_{0}′=−*J*_{1} and *Y* _{0}′=−*Y* _{1} are available. The terms of order *δ* in equation (5.3) give
5.7These expressions may be simplified by means of the Hankel approximations, valid for argument *z* that is not too small,
5.8The right-hand sides here are equal, so that in equation (5.7) we may take *Y* _{0}(*j*_{0n})≃*J*_{1}(*j*_{0n}) to obtain
5.9When convenient, the values of *K* at the critical frequency will be denoted *K*_{N}. For the parameter values in the plots, the computed values of *K*_{N} are 0.5290,1.1972,…. For comparison, approximation (5.7)_{2} gives 0.5286,1.1963,…, and approximation (5.9)_{2} gives 0.5288,1.1964,…. Hence the accuracy of the approximations is excellent.

### (b) Group velocity at the critical frequency

We now calculate the group velocity d*Ω*/d*K* at the critical frequency *Ω*=*Ω*_{N}, i.e. where *ν*=0 and *K*=*K*_{N}. As was demonstrated at the end of §2*c*, and is evident in figure 1*c*, the dispersion curve is smooth here, even though the expression for *ν* contains a square root and the value of *ν* changes from real to imaginary when *Ω* passes through *Ω*_{N}.

The most direct way of calculating d*Ω*/d*K* from a dispersion relation is to use the formula , in which subscripts denote partial derivatives. This method is impractical here, because in the dispersion relation (2.15) or (2.16), the partial derivatives of *ν* become singular at *ν*=0, and complicated cancellations and limits must be fully accounted for. Instead, it is better to write equation (2.16) in the functional form *D*(*ν*^{2},*K*/*Ω*)=0 and first calculate *d*(*K*/*Ω*)/*d*(*ν*^{2})=−*D*_{ν2}/*D*_{K/Ω}. Application of the chain rule to the algebraic relation (2.14)_{1} between (*K*/*Ω*,*ν*^{2}) and (*K*,*Ω*) gives
5.10At *ν*=0, this is
5.11evaluated at (*Ω*_{N},*K*_{N}). With the notation *J*_{0n}=∂^{n}*J*_{ν}/∂*ν*^{n}|_{ν=0} and a dash denoting derivative with respect to the argument, the series expansion of equation (2.16) in powers of *ν* gives, at *ν*=0,
5.12and
5.13Here, quantities without carets are evaluated at *K*_{N}/*Ω*_{N}, and quantities with carets are evaluated at *ϵK*_{N}/*Ω*_{N}. Expanding equations (5.12) in powers of *δ*, and rewriting the result in terms of , leads to
5.14Here use has been made of equation (5.3) for the values *K*=*K*_{N} at which *ν*=0, and of the identities *J*_{0}′=−*J*_{1}, *Y* _{0}′=−*Y* _{1}, *J*_{01}=(*π*/2)*Y* _{0} and *Y* _{01}=−(*π*/2)*J*_{0}. Hence
5.15In the last term, we have used the value of *K*_{N}/*Ω*_{N} given by equations (5.7)_{2} or (5.9)_{2}, and returned to *δ*. The coefficients here are exact. For arguments that are not too small, *Y* _{0}≃*J*_{1}, and equation (5.15) is well approximated for all *n* by
5.16On using this equation and *K*_{N}/*Ω*_{N}≃*j*_{0n} in equation (5.11), we obtain
5.17At first sight, we ought to include here the *O*(1) term from equation (5.15), for consistency with the term *j*_{0n}. However, , and so equation (5.17) includes the dominant terms. For our parameter values, it happens that . The further approximation *j*_{0n}≃(*n*−1/4)*π* in equation (5.17) gives
5.18For the parameter values in the plots, the computed values of d*Ω*/d*K* at *Ω*=*Ω*_{N} are 0.11, 0.025,…, and approximation (5.18) gives 0.11, 0.024,….

Although these calculations of wavenumbers and group velocity are somewhat detailed, and have been give in full, they are presented to show that simple final formulae may readily be obtained for every physically significant quantity in such a problem as this. These formulae are obtainable by no more than careful use of scaling arguments and series expansions.

## 6. Conclusions

This paper shows that asymptotic approximations to Bessel functions of imaginary order make possible a range of simplifications to an apparently intractable dispersion relation. The diversity of regimes analysed is wide enough that the basic identities and approximations likely to be needed in further applications of such Bessel functions are presented. Particularly useful are Dunster's identities and approximations in §3*a*,*b*, especially the main approximations (3.8) and the small-argument approximations (involving trigonometric functions of a logarithmic term) given in equations (3.10) and (3.11). The main approximations are analogous to Debye's approximations for real order; a difference is that the main approximations are everywhere wave-like, whereas Debye's approximations are wave-like only for argument greater than order. These approximations are likely to be useful in any dispersion relation containing Bessel functions of imaginary order.

## Acknowledgements

The author thanks K. R. Khusnutdinova for suggesting an asymptotic analysis of the oceanic waveguide problem treated here. Part of the work was carried out under the auspices of the Wave-Flow Interactions network in mathematics (particularly the meeting at St Andrews in June 2010), funded by the EPSRC and organized by J. Vanneste.

- Received August 1, 2012.
- Accepted August 28, 2012.

- This journal is © 2012 The Royal Society