In-out decomposition of boundary integral equations

We propose a reformulation of the boundary integral equations for the Helmholtz equation in a domain in terms of incoming and outgoing boundary waves. We obtain transfer operator descriptions which are exact and thus incorporate features such as diffraction and evanescent coupling; these effects are absent in the well known semiclassical transfer operators in the sense of Bogomolny. It has long been established that transfer operators are equivalent to the boundary integral approach within semiclassical approximation. Exact treatments have been restricted to specific boundary conditions (such as Dirichlet or Neumann). The approach we propose is independent of the boundary conditions, and in fact allows one to decouple entirely the problem of propagating waves across the interior from the problem of reflecting waves at the boundary. As an application, we show how the decomposition may be used to calculate Goos-H\"anchen shifts of ray dynamics in billiards with variable boundary conditions and for dielectric cavities.


Introduction
The transfer operator formalism proposed by Bogomolny [1] has offered a very powerful platform for the application of semiclassical methods to complex quantum and wave problems. It has long been recognised that, in the special case of cavity problems such as quantum billiards or dielectric resonators, the transfer operator is intimately connected with boundary integral methods, which provide the most effective means of exact, numerical solution. Boasman established this connection explicitly in [2] for quantum billiards, recasting the boundary integral equations as the application of a transfer operator to the wavefunction on the boundary (in the case of Dirichlet boundary conditions) or to its normal derivative (in the case of Neumann boundary conditions). Subsequent development using Fredholm theory can be found in [3,4,5,6,7], particularly on the Dirichlet case. If the boundary conditions are more general, however, it is less clear how the connection can be made exactly, although the approaches are easily related semiclassically (see [8,9,10,11], for example). In this paper we propose that a decomposition of the wavefunction at the boundary into incoming and outgoing components, defined in detail following a separation of the Green function into its regular and singular parts, provides a natural means of establishing this connection beyond leading semiclassical approximations.
Semiclassical methods originating in the field of quantum chaos have found widespread application in classical wave problems with complex or chaotic ray limits. In particular in the context of vibro-acoustics [12] or for electromagnetic wave fields [13], cavity problems with a variety of boundary conditions are prevalent. An efficient phase space flux method, the so-called Dynamical Energy Analysis (DEA), has been developed recently, predicting wave energy transport through complex structures and domains based on trajectory calculations [14,15]. A more explicit connection to boundary integral equations in terms of transfer operator methods and incoming and outgoing waves is a key ingredient for constructing hybrid DEA -wave methods [16]. This is especially valuable where semiclassically higher-order mechanisms such as diffraction or evanescent coupling are important. In addition, one often encounters the situation where semiclassical ideas provide a simple qualitative description of the problem, but where the wavelength is too large to use them for reliable quantitative predictions. A semiclassically-motivated decomposition such as proposed in this paper can then lead to an efficient implementation of fully wave-based calculations. Our treatment of the Goos-Hänchen shift is one example of this latter case; simplistic ray-dynamical models break down in the most important regions of phase space but nevertheless provide a very useful conceptual basis on which to found more accurate calculations using, for example, the method proposed here. The main technical issue to be addressed in this paper is how to decompose the wave solution at the boundary into a component approaching the boundary and a component leaving it (see Fig. 1). This decomposition is automatically achieved in the context of semiclassical approximation when the wavefunction takes the form of an eikonal ansatz, where ray direction allows us to select incoming and reflected wave components. For problems which are not semiclassical, or where the solution cannot easily be represented in eikonal form, the decomposition is less obvious. We propose that the outgoing wave ψ + (s) should be defined simply as the contribution to the boundary integral equation arising from an appropriately defined singular part of the Green function. The motivation for this definition is that, just inside Ω, the corresponding contribution to the boundary integral equation is dominated by the part of the boundary nearby and represents a wave emerging from the boundary towards the interior. This definition is consistent within semiclassical approximation with the obvious decomposition according to the direction of rays. It also reproduces expected results in the one-dimensional case and in a sense our calculation amounts to a generalisation to higher dimensions of the scattering approach normally used on quantum graphs [17,18].
There remains some ambiguity, however, arising from the possibility of defining the singular part of the Green function in different ways. The simplest definition having the correct semiclassical limit is to define the singular part to correspond to the Green function obtained by replacing the boundary locally by its tangent line or plane. The resulting decomposition can also be motivated by a transformation to momentumrepresention of the solution along the boundary. This provides a separation according to ray direction as is automatically achieved using an eikonal ansatz. We will refer to this as the primitive decomposition. The primitive decomposition is shown in [19] to allow the treatment of diffractive effects from discontinuous boundary conditions, for example, or to provide a basis for the semiclassical coupling of transfer operators between multiple domains with different local wave number.
The primitive decomposition displays singularities in momentum representation near the case of critical reflection where waves arrive at the boundary almost tangentially. Although one might still define an associated transfer operator before semiclassical approximation, such singularities will at least lead to issues such as slow convergence. The regime of near tangential angle of incidence is particularly important when describing leakage from dielectric resonators. Here the "critical line" in phase space, along which refractive escape switches on, plays an important role in emission patterns and decay rates. We therefore propose a second separation of the Green function into singular and regular parts based on a splitting of the boundary integral into distinct contours in the complex plane, one of which captures the Green function's singularity. We call this second decomposition the regularised decomposition. It successfully accounts for the important effects of boundary curvature at critical reflection and remains well-behaved when waves arrive nearly tangentially. We examine the regularised decomposition in detail for the special case of circular cavities.
We find in general that the boundary integral equations can be cast in the form where the operatorŜ depends on the geometry of the domain but is completely independent of the boundary conditions: in this paper we will refer toŜ as the shift operator. This equation is naturally interpreted as saying that ψ − is obtained by propagating the outgoing wave ψ + across the interior until it returns to the boundary as an incoming wave. For the case of circular cavities we provide explicit, closedform expressions forŜ. We emphasise that this is achieved even in the presence of variable boundary conditions such as those treated in [20,21], where the problem as a whole is nonintegrable. We show that, closing the system by using the boundary conditions to relate the outgoing to the incoming component by a reflection operator r, we arrive at an overall transfer operator taking the form of a stroboscopic map T =rŜ. Semiclassical treatments provide a simple interpretation ofr as acting on rays by a Goos-Hänchen shift, so that the overall operator offers a simple explanation of the Goos-Hänchen-perturbed ray dynamics of the type used in [22,23,24,25] to treat dielectric cavities. Furthermore, the regularised decomposition allows a more complete treatment of the case of critical reflection (or refraction), albeit in that case with a less simple interpretation in terms of rays. At a formal level, the construction of a transfer operator as a map between boundary functions is similar to the scattering-matrix approach developed by Smilansky and coworkers [26,27,28] or the transfer-operator approach by Prosen [29,30,31]. Similar to our approach, a wave solution is decomposed in those calculations into components passing in either direction through a surface of section cutting through the domain of interest or split into an incoming and outgoing wave component at (parts of) the boundary of a cavity. The formalism by Smilansky et al and Prosen rely on being able to construct explicitly a scattering matrix for the regions lying to either side of a surface of section. They are thus most natural when decomposing the problem in hand into two 'half-problems' separated by the surface of section. Such a treatment becomes less obvious when using as the section the full boundary of a closed cavity, the most natural surface of section when starting from boundary integral equations. It is worth mentioning that the transfer operator construction proposed by Prosen [29,30,31] is similar to the 'primitive decomposition' introduced above. Prosen developed a more general approach dealing also with smooth potentials, but the restricted circumstances we assume here allows a simpler derivation and a generalisation to a regularised decomposition which deals better with the case of critical reflection, for example.
The emphasis in this paper is quite different, however. We argue that, starting from the boundary integral equation, one can naturally decompose the Green functions occurring there into singular and regular components that are closely related to the incoming and outgoing component of the boundary wave function. The regularised decomposition provides in particular a description which stays regular at grazing incident angles and which makes it possible to distinguish propagating and evanescent (exponentially decaying) wave contributions without explicitly referring to a (scattering) channel basis.
The regularised decomposition approach applies in it simplest form to cavities that have convex, analytic boundaries. The boundary needs to be analytic so that the contour integration methods used to define singular and regular parts of the Green operators can be defined. Assuming that the boundary is convex allows a more direct leading semiclassical description of the shift operatorŜ in terms of the Green operators, without the need to account explicitly for the cancellation of ghost orbits (see [7]). Treatments of piecewise analytic boundaries, with a more explicit account of corners, and of nonconvex billiards will be the subject of future investigation.
We conclude this section by summarising the content of the paper. We begin in Sec. 2 by establishing notation and outlining the main technical features of the calculation. The use of the decomposition of Green operators into local and nonlocal parts to define incoming and outgoing waves is illustrated in Sec. 3 for the onedimensional case, where calculations are elementary. This is generalised to twodimensional cavities in Sec. 4 using the primitive decomposition. The regularised decomposition, taking more explicit account of boundary curvature, is treated in Sec. 5. Finally some model two-dimensional cases are examined in Secs. 6 and 7, providing in particular explicit analytic descriptions of the shift operators for the circle, and conclusions are offered in Sec. 8.

Overview
We now summarise the main features of our calculation, concentrating on the case of two-dimensional cavities to fix notation. We emphasise that the main ideas in this section generalise to higher dimensions.
Let ψ(x) satisfy the Helmholtz equation in the (two-dimensional) domain Ω. We use the same symbol ψ(s) to represent the restriction of the solution to the boundary ∂Ω, on which s is an arc-length coordinate, and denote the normal derivative by µ(s) = ∂ψ ∂n .
. The boundary integral equation for ψ may then be written formally as where we denote the Green operatorŝ in which x denotes a generic point in the interior of Ω and s and s label points on the boundary. Ω s ψ ψ x + − Figure 2. The wave emerging from the part of the boundary near x, as x → x(s), provides a natural definition of the outgoing component ψ + . We define this more directly as the contribution from the singular part of the Green function. The regular remainder defining ψ − represents the contributions originating from the rest of the boundary, arriving at s having passed through the interior.
The key step in the transformation we propose is to perform a basis change in which we replace the boundary wavefunction ψ(s) and normal derivative µ(s) by incoming and outgoing components ψ − (s) and ψ + (s), defined in the simplest implementation so that Motivated by the schematic picture of Fig. 2, we argue that ψ + (s) is naturally defined as corresponding to the part of the boundary integral in which s approaches s. This should capture in particular the singularity of the 2D Green function, as s → s. We therefore separate the Green operators into "singular" and "regular" parts,Ĝ (with an analogous decomposition ofĜ 1 ), whereĜ sing 0 is an operator which should give the contribution to the boundary integral from s near s and should in particular account for the singularity (6). We then define and Note that, in the schematic picture of Fig. 2, ψ − (s) now collects the contributions to the boundary integral equation from the remainder of the boundary, arriving at s having crossed the interior of the cavity. This is consistent with our requirement that ψ − should represent the part of the solution that is incoming at the boundary. We discuss in the coming sections the detailed calculation that arises once the decomposition (7) has been defined more explicitly. Here we summarise the main conclusion, which is, that the boundary integral equation (2), expressed in terms of ψ + and ψ − , can be recast formally as equation (1). We emphasise that this equation derives entirely from the boundary integral equation and has nothing to do with the boundary conditions. For a given domain Ω, we get the same shift operatorŜ in (1) regardless of whether Dirichlet, Neumann or any other boundary condition is imposed.
The boundary conditions provide us with a second, independent relationship between ψ − and ψ + . Suppose that the boundary conditions are prescribed in the form of a linear relationship between ψ and µ on ∂Ω. (This assumption must be modified for the treatment of dielectric and other problems where coupling to an exterior solution is involved, as described in Sec. 7). We re-express this in terms of the functions ψ − (s) and ψ + (s) as where we callr the reflection operator, for obvious reasons. Combined with (1), we find is the transfer operator for the problem as a whole. Eigenvalues are found by solving the secular equation We find indeed that, when semiclassical approximations are used forŜ andr, we recover the usual transfer operator expected from [1] (up to a simple conjugation of the operators involved: see Appendix B). Note that the transfer operatorT is presented as a stroboscopic map, composing the shift operatorŜ, determined entirely by the shape of Ω, with the reflection operatorr, accounting for boundary conditions. In semiclassical approximation, S is determined by the usual ray dynamics assuming specular reflection at each bounce. The reflection operator can be accounted for in this picture simply by adding appropriate phases as rays are reflected from the boundary. An alternative point of view, which has been found to be fruitful in the treatment of dielectric problems, for example, is to regard the variable reflection phases inherent inr as adding a perturbation to the usual ray dynamics in the form of a Goos-Hänchen shift [22,23,24,25,32,33]. In our approach, this is easily understood from the stroboscopic nature of the full operatorT . Writingr formally as an exponential and using an analogy with quantum mechanics in which k plays the role of 1/ , we can interpretĥ as the generator ofr. The corresponding effect on rays is obtained simply by replacingĥ by a symbol h(s, p) on the boundary phase space and using this as a classical Hamiltonian. Note that, because the reflection phases are formally of O(1) in an eikonal expansion in 1/k, the Hamiltonian h(s, p) so defined is small and therefore only perturbs the main ray dynamics (and that is why we are equally entitled to leave the ray dynamics alone and incorporate reflection phases instead in the next-to-leading, amplitude-transport phase of the eikonal expansion [34]).

The one-dimensional case and quantum graphs
In order to motivate equations (8) and (9) further as a means of decomposing the boundary solution into incoming and outgoing components, we consider now the onedimensional case. The most direct analogue leads very simply to the standard picture one uses for the wave dynamics on quantum graphs [17].
one finds that the Green operators can be represented by the 2 × 2 matrices and where I is the identity matrix and will be found to play the role of the shift operator.
In this one-dimensional case, instead of separating the Green function into a singular and a regular part, as in (7), we separate the diagonal from the nondiagonal part That is, the role played in higher dimensions by the singular part (as s → s) of the Green function is in one dimension taken over by the diagonal part (in which s = s). Then, defining and in analogy with (8) and (9), we find that the boundary integral equation (2) becomes Note that in the case of quantum graphs we simply apply the same calculation independently to each bond and find an overall shift operator with diagonal blocks of the form (15). The problem is closed by applying boundary conditions at vertices to provide scattering operators that generaliser in (10) and couple the bond-labelled blocks ofŜ.

The primitive decomposition
The simplest decomposition in the two-dimensional case is obtained by defining the singular part of the Green operator to be the result of replacing the boundary locally by its tangent line, as illustrated in Fig. 4.  Before taking the limit, the kernel of the full Green operator entering (3) takes in 2D the form where L(x, s ) denotes the distance between the generic interior point x and the boundary point s and H 0 (z) is the outgoing Hankel function. If the boundary is replaced by the tangent line at s, (see Fig. 4), this distance can be written where ε denotes the normal distance from x to s. We define a corresponding singular component of the Green operator bŷ The boundary function µ(s) entering the integrand in (19) and in analogous equations below is here understood as being the periodic extension of the true boundary function onto the real line with period equal to the perimeter of ∂Ω. (This is understood to hold both for open and closed boundaries in order to unify notation.) This operator takes a remarkably simple form when expressed in a Fourier basis (or, in quantum mechanical language, a momentum representation): formally defines a normal momentum operator and p = 1 ik ∂ ∂s denotes a tangential momentum operator. More precisely, the action of the operator P on a function ψ(s) is defined with respect to the Fourier transform to beP where integration is over a contour in the complex p-plane that passes above the branch singularity p = −1 and below the branch singularity at p = 1. The singular part of the operatorĜ 1 is even simpler: whereÎ is the identity operator. These relations are derived explicitly in Appendix A and reconfirm the well-known jump condition important to boundary integral methods more generally. The outgoing wave defined in (8) can now be written providing a two-dimensional generalisation of (17). Together with the condition , which is what one might expect intuitively on the basis of a Fourier representation of the solution near the boundary: a boundary solution ψ(s) = e ikps extends to an interior solution ψ ± (x) exp[ik(ps ∓ 1 − p 2 n)] near the boundary, where n represents a normal coordinate increasing towards the interior and the sign is determined by whether the extended wave moves towards or away from the boundary. Taking a derivative with respect to n then gives the form for µ written above.
The relevant semiclassical expressions are described more explicitly in Appendix B. Note also that, although (28) requires an operator inversion, the operator to be inverted is close to the identity in semiclassical approximation for convex cavities. One can then easily incorporate the inversion into higher-order approximations in that case.

Regularised decomposition for analytic boundaries
The primitive decomposition leads to singularities in momentum representation around p = ±1 that are unphysical for curved boundaries (see the discussion in Sec. 6.2 around Fig. 7, for example). This singular behaviour can be removed by taking account of boundary curvature more explicitly. We describe in this section an alternative decomposition which achieves this goal (in two dimensions). As a concrete example, we will treat circular cavities in detail in Sec. 6, for which all the relevant operators can be fully understood analytically. We emphasise, however, that much of the discussion generalises to other geometries with analytic boundaries.

Decomposition by integration contour
To motivate the construction, let us examine the action of the Green operatorĜ 0 , on the boundary plane wave In asymptotic analysis it is natural to treat the resulting integral in which denotes the length of the perimeter, using the method of stationary phase and, where appropriate, exploiting the asymptotic representation of the Hankel function, which begins so that 2πikL(s, s ) .
One can identify in particular two kinds of contribution: those from the endpoints at s = s and s = s + , and those from saddle points. At leading order, the asymptotic contribution of endpoints is obtained by Taylor-expanding the chord-length function and truncating at the linear term, so that we use respectively for the two endpoints. We now point out that this truncated length function is precisely the length function (following the limit ε → 0) of the tangent line used in the previous section to defineĜ sing 0 . That is, the primitive definition ofĜ sing 0 µ p (s) amounts essentially to a separating out of the leading-order endpoint contributions from the integral in (30). This suggests that a better definition ofĜ sing 0 might be obtained by including higher-order corrections to the leading endpoint contribution already used. This would account in particular for curvature, which enters the Taylor series for L(s, s ) at third order. We assume in this section that the boundary is analytic. Then a systematic treatment of the full asymptotic expansion originating with the endpoints is obtained by promoting the integral in (30) to the status of a contour integral in the complex s -plane and then isolating the contour component Γ 0 emerging from the endpoints and following a path of steepest descents. For appropriate test functions we might then define, simply, This definition has the advantage of automatically offering the systematic asymptotic expansion of the endpoint contribution while also permitting, in principle, exact analysis. In practice, this definition is complicated by the fact that the topology of Γ 0 depends in general on the test function. In particular, for the case of plane waves µ p (s), Γ 0 may take one of two topologically distinct routes, depending on whether p 2 < 1 or p 2 > 1 as illustrated in Fig. 5 for the case of the circle. The contour for propagating waves with p 2 < 1 separates into a two-part endpoint component Γ 0 -one part starting at s = 0 and ascending into the upper-half plane and the other part ascending from the lower-half plane and ending at s = s + -and a separate component Γ 1 which completes the contour and caters for real saddle points corresponding to the rays passing from s to s across the interior. Note that such propagating waves correspond in the language of the scattering approach [26,27,28,29,30,31] to open modes, and to reflect this we let H o denote the subspace of boundary functions defined by p 2 < 1 in momentum representation.
Thus, we defineĜ sing 0 on H o by (31) and by fixing the contour Γ 0 to be the lines ascending vertically from s = s to the upper-half plane and ascending vertically from the lower-half plane to s = s + . Note that by fixing Γ 0 to rise vertically it will in general deviate from the strict steepest-descent contour but will give the same result up to the exponentially small differences. These may arise when the respective contours go around singularities on different sides. Using a fixed contour has the advantage of letting us define the actionĜ sing 0 on H o without reference to a specific basis within that subspace.
IfĜ sing 0 is defined by (31), the action ofĜ reg 0 on H o is then defined bŷ where Γ 1 is now also fixed and understood to descend vertically from the upper-half plane to s = s, to move across the real axis to s = s + and then descend vertically into the lower half plane. Note that Γ 1 can be deformed so that it avoids entirely the endpoints of the complete integral and is dominated in asymptotic analysis by saddle points corresponding to interior-crossing rays (and potentially their complex generalisations). We reiterate that at leading order in asymptotic analysis these definitions reproduce the action of the primitive decomposition defined in the previous section.
For the subspace H c spanned by closed modes with p 2 > 1, asymptotic analysis uses a different contour decomposition, reflecting the movement of contributing saddle points into the complex plane, directly below s = s (for p < −1) or above s = s + (for p > 1). In this case, the decomposition of contours used to defineĜ sing 0 andĜ reg 0 is not as obviously defined, suggesting that there are alternative, a priori equally sensible decompositions of closed boundary modes into incoming and outgoing components. In fact, the asymptotic form of the integral (31) is in the middle of a Stokes transition when p 2 > 1 is real. Illustrated in Fig. 5(b) is the contour decomposition for a test function with p < −1 in the circle. Here the endpoint contour segment Γ 0 descends vertically from s = s and runs into a saddle. If p moves to either side of the real axis, Γ 0 passes to one side or other of the saddle. We may then either include the saddle separately, by including a contour component Γ 1 running over it completely (see Fig.  C1 in Appendix C) or miss it entirely. We choose the latter option, which amounts to choosing Γ 1 = 0 and letting Γ 0 account for the entire contour, as in Fig. 5(b). This implies the choiceĜ reg 0 H c = 0. See Appendix C for discussion of the alternative convention.
The Green operatorĜ 1 could be decomposed analogously by separating the appropriate boundary integral using the same contour components. We note that, when taking the limit x → s in this approach, a jump condition analogous to (24) is satisfied, but that there remains a nontrivial correction: whereK is obtained by putting x directly onto the boundary and α is an incidence angle defined in Fig. B1. So far our discussion has suggested that each ofĜ 0 andĜ 1 might be decomposed individually by separation of contours. This approach has the advantage of being conceptually simple while allowing systematic asymptotic approximation beyond the leading primitive terms, and also allowing exact analysis in principle. As will be shown in Sec. 5.2 below, it may be advantageous in practice to adopt a variation of this approach. In particular, the calculation of the shift operatorŜ is greatly simplified if the identityĜ is satisfied. Note that this identity holds for convex domains at leading order for any of the decompositions we consider in this paper (including the primitive decomposition of section 4). It is also shown in section 6.1 that the definitions ofĜ sing 0 andĜ sing 1 so far given yield the identity (35) exactly in the case of the circle.
More generally, an alternative approach is to guarantee that (35) holds by choosing one of the following: A let (31)  Because (35) holds in any case at leading order (for convex domains), higher corrections can be systematically included in an asymptotic analysis of this alternative approach. In this context, option B in particular is attractive becauseĜ sing 1 has the simple leading approximationĜ sing 1 −Î/2. Finally, we note that if the length function L(s, s ) suffers singularities in the complex plane (in addition to those at s = s and s = s + ), the functionĜ sing 0 µ(s) defined by (31) may exhibit discontinuities whenever these singularities cross Γ 0 (as s is varied). Any such discontinuities will be small, however, as the integrand decays exponentially along Γ 0 . In practical terms, such singularities can be circumvented by not allowing the contour segments to extend indefinitely into the complex plane; that is, we may truncate Γ 0 (and therefore Γ 1 ) at a finite distance from the real axis, so that passing singularities are missed. Note that the detailed decomposition will then depend on the heights at which the contours are truncated, but only to the extent of exponentially small corrections in asymptotic analysis. For the case of circular cavities described more fully in section 6, there are no such difficulties and we can give a complete analytical description of the resulting decomposition.

Shift operator for the regularised decomposition
Let us now outline how a decomposition of the Green operatorsĜ 0 andĜ 1 into singular and regular parts leads to a restatement of the boundary integral equations in terms of a shift operator. We do this initially without making any particular assumptions about the detailed properties ofĜ sing 0,1 andĜ reg 0,1 and then point to the simplified results that may be obtained if explicit properties such as (35) are assumed.
We start, in general terms, with equations (8) and (9) defining the outgoing and incoming components of the boundary solution. From (8), we may eliminate µ = Ĝ sing 0 −1 ψ + +Ĝ reg 1 ψ and, using also ψ = ψ + + ψ − , we then find from (9) that Note that the primitive version (27) is obtained as a special case of this equation.
If we now also assume that the singular componentsĜ sing In other words, the shift operator is then simplŷ In fact there is one further, very suggestive simplification in this case. Consider the exterior problem, in which boundary data (ϕ(s), ν(s) = ∂ϕ/∂n) are once again supplied on ∂Ω but where the Helmholtz equation is now assumed to apply in the exterior region Ω , with radiating boundary conditions imposed at infinity. We change the notation for this boundary data to emphasise that the solution of this exterior problem is different to that of the interior one. Then, assuming that the normal is still defined to be pointing out of Ω (and into Ω ) the boundary integral equation replaces (2). HereĜ 1 denotes the same operator used elsewhere in this paper -in which the limit x → s in (4) is taken from the interior. The main differences from (2) are some sign changes to reflect that the normal points into Ω and the addition of an identity operator to the last term to correct the jump condition. DefiningĜ sing 0,1 andĜ reg 0,1 exactly as in the interior case, and letting the decomposition ϕ = ϕ − + ϕ + be such that we get the analogue of (36). In this case, if we assume thatĜ sing 0 andĜ sing 1 are chosen so that (35) holds, we find that In other words, imposing (35) means that the regularised decomposition correctly identifies the component on the boundary that extends to the outgoing component in the farfield.

Further notation for the regularised decomposition
We now suggest notation that makes a more direct analogy with the notation already used for the one-dimensional problem and for the primitive decomposition. This notation also simpifies the representation of boundary conditions as reflection operators in the model problems of Secs. 6 and 7. Let us define operatorsD 0 andD 1 bŷ That is,D 0 andD 1 are simply scaled and shifted representations of the singular parts of the Green operators. We further define operatorsŜ 0 andŜ 1 bŷ These equations replace (13) and (14) for the one-dimensional calculation and (25) and (26) for the primitive decomposition. In particular, the primitive version is retrieved by replacingD 0 → 1/P andD 1 → 0. Note also that we can alternatively writê which may be distinct if (35) is not imposed, as is the case for the primitive decomposition. We confine our attention once again to the interior problem. Then the outgoing and incoming components of the boundary solution are respectively and Eliminating ψ and µ in favour of ψ − and ψ + in these equations leads to (36), in which the shift operator takes the form This generalises the primitive equation (28). If identity (35) is imposed then this reduces to the much simpler form previously provided in (37). These operators are evaluated explicitly for the case of a circular domain in the next section.

Applications to a model problem -circular cavities
In order to illustrate the application of the in-out decomposition to a nontrivial problem, we consider now the case of the circle with variable Robin boundary conditions of the form with This includes as a special case a = 0 the examples treated in [20,21]. Such problems provide an interesting challenge for our approach when regions of the boundary exist for which F (s) > 0. This is because the system then supports modes that are localised in the full 2D problem near the boundary [20,21] and that have on the boundary itself significant components in momentum representation with p 2 > 1. The regions around and beyond critical lines p 2 = 1 then play a significant role and allow us to test decompositions across the transition from open to closed modes. We emphasise that, although the Green function decomposition can be achieved analytically in this geometry, the variable boundary conditions make the complete problem nonintegrable [21].

Shift operator for the circle
We begin by finding the shift operator for the circle. In this case we can write explicit analytical results for the various operators defined in Secs. 5.2 and 5.3. Furthermore, we find that identity (35) automatically holds whenĜ sing 0 andĜ sing 1 are respectively defined by (31) and (33).
The chord-length function for a circle of radius R is Then an application of Graf's theorem [35] allows us to write the kernel of (3) in the form where z = kR. In a Fourier basis (or momentum representation),Ĝ 0 is then represented by a diagonal matrix whose entries are The two alternative forms given here for (Ĝ 0 ) mm respectively provide the decomposition ofĜ 0 for |m| > z and for |m| < z. That is, one finds (although we do not show the detailed integrations here) following the recipe described in Sec. 5.1 that the singular and regular parts of the Green operatorĜ 0 are represented by diagonal matrices with entries Similarly, one finds that has the decomposition We also find that In summary, we have shown that (S) mm = (S 0 ) mm = (S 1 ) mm holds in the circular case (as does (35)). It is also useful to record that the scaled singular parts defined in Sec. 5.3 have diagonal elements and primes indicate derivatives with respect to z. We note that the leading Debye approximation for H m (z) [35] yields, for z 1, consistent with the primitive decomposition on making the identification p = m/z, while is asymptotically of higher order as z = kR → ∞. Such primitively defined operators have been successfully used in [19] to describe eigenfunctions of a disk with boundary conditions changing discontinuously from Dirichlet to Neumann, where diffractive effects are important. We show in the more detailed calculation of Sec. 6.2, however, that they lead to an in-out decomposition with undesirable features in the evanescent region of momentum space (see Fig. 7).

Circular cavity with variable Robin boundary conditions
Having characterised the shift operator in the previous section, the remaining step is to formulate the prescribed boundary conditions in terms of a reflection operator mapping the incoming to the outgoing boundary solution. To this end we simply substitute the boundary condition (45) into the equation (43) that defines the outgoing wave ψ + , which leads to Î − iD 1 − iD 0F ψ + = Î + iD 1 + iD 0F ψ − . Hence whereX =D 0F +D 1 .
Note thatX is similar to the operator X =D the evaluation is straightforward here, requiring simply the determinant of a tridiagonal matrix for the special case (46). Alternatively, the secular equation can be written directly as a difference equation as described in [20,21], albeit expressed in a different basis to that used here. We note that the form of the reflection operator in (48), is independent of the choice ofF in (45). Our decomposition in terms of shift and reflection operators thus allows us to treat general boundary conditions easily once the shift operator has been constructed. (Of course, the transfer operator does not have tridiagonal form in general.)

Evanescence in solutions
We now describe some particular solutions of this model problem, concentrating on their behaviour in the region of closed modes p 2 > 1.   [20,21] in terms of waves localised near the boundary wherever F (s) > 0. We now offer a brief summary of this explanation, borrowing from the discussions in [20,21].
Consider first the case where F > 0 is a constant and kR 1. Then one can find quasimodes decaying exponentially into the interior (denoting p = m/kR and s = Rθ) of the form ψ(r, θ) e −kq(R−r)+imθ = e −kq(R−r)+ikps , satisfying the interior wave equation if p 2 = 1 + q 2 and the boundary conditions if q = F . Note that such boundary-localised solutions occupy the closed subspace corresponding to p 2 > 1. If F now varies, but over a scale longer than that of a wavelength, we should therefore expect localised solutions formed by deformations of these states, spread in momentum representation over the region 1 + F 2 min < p 2 < 1 + F 2 max , where F max = max(max(F (s)), 0) F min = max(min(F (s)), 0). Exact eigenfunctions will then generically couple to such local solutions and lead to the plateau structures shown in Fig. 6. Heuristically we can interpret the observed plateaux as being the result of tunnelling in momentum space from the primary quasimode around m = 33 to an "allowed region" 1+F 2 min < p 2 < 1+F 2 max introduced by the boundary-localised mode. In fact, the four examples in Fig. 6 were chosen to have the same values of F min = 0 and F max = 1, and one does indeed then observe that the plateaux occupy the same interval, indicated by the shaded strip in the figure. . The outgoing boundary intensity |ψ reg + | 2 defined by the regularised decomposition is shown as a heavy continuous curve over the plateau region for case (c) of Fig. 6: as in Fig. 6, the "allowed region" 1 < p 2 < F 2 max is indicated by the shaded strip. Also shown, as a thin continuous curve, is the version |ψ prim + | 2 obtained from the primitive decomposition. The thin dashed curve shows the primitive incoming intensity |ψ prim − | 2 . The primitive and regularised outgoing intensities differ significantly only near m = kR ≈ 74, where the primitive case shows an expected divergence. The smaller, but significant, incoming intensity observed in the primitive case is qualitatively different from the regularised case, where the incoming intensity is zero or exponentially smaller than the outgoing intensity (depending on the conventions used to define the contour Γ 0 in Sec. 5.1).
This example offers a useful illustration of the most significant differences between the regularised and primitive decompositions -and their corresponding definitions of ψ ± (s). In Fig. 7 we have selected the eigenfunction corresponding to case (c) of Fig. 6 and shown it (as the thick continuous curve) in more detail over its plateau region 1 < p 2 < F 2 max . Also shown are the primitive version (as a thin continuous curve) and the primitive incoming solution (as the thin dashed curve). As expected, the regularised and primitive outgoing intensities are very close, hardly distinguishable on this log plot, except that the primitive intensity becomes significantly larger in the region of critical incidence p 2 ≈ 1. We emphasise that this divergence is an artefact of the primitive decomposition, and its failure to deal properly with curvature at critical reflection, rather than an intrinsic characteristic of the solution itself.
The second difference is that, in the case of primitive decomposition, we obtain a smaller, but nonetheless significant, incoming component even in the "evanescent" region p 2 > 1, see Fig. 7. By contrast, the incoming component vanishes for the regularised decomposition. The nonvanishing primitively-defined incoming wave can be understood further by explicitly evaluating the shift operator for the primitive decomposition. Evaluating expression (28) for the special case of the circle yields a primitive shift operator represented by a diagonal matrix with entries which needs to be compared withŜ in (47) for the regularised decomposition. For closed components with m 2 > z 2 ≡ (kR) 2 , and using the Debye approximation for the numerator, Eq. (50) vanishes at leading order. However, the next-to-leading order contributions do not cancel. Therefore ψ prim − =Ŝ prim ψ prim + does not vanish as does the regularly-defined incoming wave, or even decay exponentially as one would expect for properly-defined evanescent waves. Instead, it is only smaller than ψ prim + by a factor scaling algebraically with 1/k. This failure of the primitive decomposition to separate evanescently incoming and outgoing solutions at all orders puts it at a significant disadvantage at describing such solutions compared to the regularised decomposition.

Husimi functions and Goos-Hänchen shifts
A particularly helpful means of presenting eigenmodes of this problem is by the boundary Husimi functions. Here we are guided by previous applications to dielectric problems [36,37,38,39], where a definition of boundary Husimi functions has been proposed. This amounts in our language to applying the conventional definition of Husimi functions to the primitive in-out components discussed further in Appendix B. In this paper we choose to work instead with a variation in whichP 1/2 is replaced by its regularised analogueD where the state |s 0 , p 0 corresponds to the boundary function and the inner product is Using (52) rather than (51) gives a Husimi function that is better behaved around the critical line p 2 = 1 but is otherwise very similar in the region p 2 < 1 and does not have a significant impact on the discussion following. A particular pair of Husimi functions H + (s, p) and H − (s, p) calculated in this way is shown in Fig. 8, corresponding to the eigenmode labelled (d) in Fig. 6. The incoming and outgoing modes in the figure are qualitatively similar but display differences in detail that are explained by the application of the reflection operator in (48) -or equivalently by the Goos-Hänchen (GH) shift [22,23,24,25,32,33]. Goos-Hänchen and related perturbations of ray dynamics are typically explained in terms of a shifting of the centre of a Gaussian or other localised beam following reflection from a boundary or interface. Such approaches provide a clear physical interpretation of the effect but seem unnatural in the context of calculating eigenfunctions, where no Gaussian beams appear in the rather more complicated wave patterns that one finds in typical cases. We now point out that the in-out decomposition of boundary solutions explored in this paper provide a much more direct way of describing the effect in this context. It makes it possible to explain the GH shift in terms of the tools that are naturally used to calculate eigenfunctions with chaotic or otherwise nonintegrable ray dynamics. Furthermore, the improved treatment of critical reflection such as provided by the regularised decomposition may allow improved descriptions of effects such as Fresnel filtering that are observed in associated parts of phase space, although we do not pursue that issue explicitly in this paper.
The incoming and outgoing solutions respectively satisfy the equationŝ whereŜ describes the ray dynamics of specular reflection and is common to all cavity problems occupying the same domain, irrespective of boundary conditions. The GH shift may be incorporated simply by representingr in the exponential form (11) and regarding it as a quantisation of the flow of a Hamiltonian generator h(s, p). The complete dynamics appropriate to ψ ± are then obtained simply as a stroboscopic intertwining of this "Goos-Hänchen map" and the regular, specular ray dynamicsalbeit in opposite orders for ψ + and ψ−.

The reflection operator in (48) is generated bŷ
Away from the critical line p 2 = 1 this has the leading-order "primitive" symbol 1 − p 2 − iF (s) and the Goos-Hänchen map can be approximated at leading order by Note that ∆s and ∆p are O(1/k) and therefore act as perturbations of the usual specular ray dynamics. Note also that there is a nontrivial shift in momentum as well as position here because the local reflection phase depends on s. In the case of planar reflection, where the reflection phase is constant along the interface, the Goos-Hänchen map changes s but not p, as prescribed in classical descriptions of the effect. The breaking by the boundary conditions of integrability is strong enough in the example of Fig. 8 to localise the mode within an island chain of Goos-Hänchenperturbed ray dynamics (see [23] for a treatment of this effect in the context of dielectric cavities). In Fig. 9 we illustrate the Goos-Hänchen vector field (∆s, ∆p) over part of phase space surrounding the point (s, p) = (0, 1/2), which lies at the centre of the middle-upper lobes of the Husimi functions of Fig. 8: this is where the incoming and outgoing Husimi functions differ most in this example. One sees that the vector field is indeed entirely consistent with the deformation of the incoming Husimi function relative to its outgoing equivalent: features lying at the bases of the Hamiltonian vectors for the former case lie at their tips in the latter case. Note that Poincaré plots of the corresponding dynamics (Goos-Hänchen map following the specular map in the outgoing case and preceding it in the incoming case) show a similar deformation (not shown).

Application to dielectric cavities
Our final application is to the problem of calculating resonant modes in dielectric cavities. The problem is solved explicitly here for the circular cavity with boundary conditions appropriate to TM polarisation. This simple problem can be solved analytically by elementary means: our aim, however, is to illustrate the application of the method to inside/outside problems in a way that it might later be generalised to less trivial, noncircular domains.
In this case, we must separately solve two wave problems; the interior one whose solution and normal derivative we denote as usual by ψ and µ respectively, and which is such that x ∈ Ω, and the exterior problem whose solution and normal derivative are respectively denoted by ϕ and ν and for which Here the refractive index inside the cavity is n > 1 while outside we set n = 1.
Following the discussion of Sec. 5.2, the interior and exterior problems separately define decompositions ψ = ψ − + ψ + and ϕ = ϕ − + ϕ − and shift operators such that The shift operatorsŜ int andŜ ext differ by using Green functions corresponding to different values of the refractive index and by the boundary integral equation for the exterior problem having a unit normal pointing into rather than out of the domain concerned, (see Eqn. (38)). According to the conventions suggested in Sec. 5.2, where the endpoint contour Γ 0 is chosen for closed modes to account for all of the boundary integral, we find that the external shift operator vanishes, S ext = 0 and Eqn. (40) holds (see Appendix C for discussion of alternative conventions). Of course it is the field ϕ(x) in the exterior x ∈ Ω that is physically observed here: see Appendix D for a discussion of the determination of the full solution, inside and outside the cavity, from its restriction to the boundary. The end result will be presented more compactly if we use the shift operator for the external problem to define a Dirichlet-to-Neumann mapF , such that In general (after some manipulation), and we note thatF itself is independent of conventions used in the in-out decomposition. Note also thatF is determined by the condition of having no incoming waves at infinity and is independent of boundary conditions on the cavity itself. We next apply the cavity's boundary conditions. To be concrete, we impose TM boundary conditions ψ(s) = ϕ(s) and µ(s) = ν(s), but emphasise that the following discussion generalises easily. Together with (54), this provides a boundary condition formally very much like that of Sec. 6 and leads to a similar-looking reflection operator except that, in this case, . Note that this reflection operator can alternatively be written This allows for a more direct comparison with classical expressions for the corresponding Fresnel reflection coefficient r = n cos α − cos β n cos α + cos β .
This comparison is achieved on making the identifications ∼ cos β ≡ i n 2 p 2 − 1 appropriate to the ray-dynamical limit, where α and β are respectively the angles of incidence and refraction of a ray hitting the boundary from inside.
We can now determine resonant modes as solutions of the secular condition This is once again formally much as in the model example of Sec. 6. The main qualitative difference is that, whereas the restriction of the reflection operator to the open subspace has eigenvalues on the unit circle in the example of Sec. 6, here the relevant eigenvalues ofr lie inside the unit circle. reflecting the lossiness of the system (to radiation to infinity). As a consequence, the solutions of the secular equation are obtained here for complex values of k, as should be expected for resonant modes with a finite lifetime. We point out that composing the shift operatorŜ with a subunitary reflection operatorr as in (57) provides a very natural wave analogue of the ray-dynamical approach taken, for example, in [40,41,42,43]. There, emission from dielectric cavities is treated by iterating initial densities of rays under a composition of the raydynamical analogues ofŜ andr, namely a Perron-Frobenius operator encoding the Poincaré-Birkhoff map of the internal dynamics and a transmission loss determined from classical Fresnel coefficients. The combined operator in (57) simply replaces each of those elements by a corresponding wave map. It therefore offers a basis for treating inherently wave effects within the broader approach of [40,41,42,43], particularly, for example, the more complicated transmission losses encountered around the critical lines (see also further discussion below).
We now specialise these results to the case of a circular cavity. One then finds from (55) thatF is diagonal and has elements (Of course this can also be deduced directly from (54) by elementary means but our aim here is to offer a discussion that generalises.) We also find that (for all m) and from which the reflection operator follows by (56). Note that this provides a means of calculating curvature corrections to Fresnel laws, such as described in [44], while also incorporating the phase information relevant to, for example, Goos-Hänchen shifts. The elements (F ) mm undergo a transition across the critical lines of the exterior problem, defined by |m| = kR. They are approximately imaginary for |m| < kR and define through (56) diagonal elements of the reflection operator such that |(r) mm | < 1: corresponding resonant modes leak directly by refraction to the exterior and are shortlived. Whispering gallery modes lie within the open subspace of the interior problem, for which |m| < nkR, and the closed subspace of the exterior problem, for which |m| > kR. The corresponding elements ofF are approximately real but retain small positive imaginary parts and lead to diagonal elements ofr which are close to, but slightly less than, unity in magnitude: the corresponding whispering gallery modes therefore have a long but finite lifetime.
A representation ofr in exponential form (11) defines a generatorĥ with significant nonHermitian components in the region |m| < kR and smaller but nonvanishing nonHermitian components over kR < |m| < nkR, reflecting the lossiness of the system as a whole. Away from the critical lines |m| = kR, where the reflection phase is not rapidly varying, we may approximate this operator straightforwardly by a Goos-Hänchen map for which the Hamiltonian generator h(s, p) has imaginary components and is obtained as the ray-dynamical symbol ofĥ = i log(r)/k. Across the critical lines |m| = kR, where the reflection phase is rapidly varying, the raydynamical shifts corresponding tor are less easy to work out. However, one can still formally user to define a ray-dynamical map -it is just not as close to the identity map and its generator is not as easily approximated. Further investigation of this issue will hopefully lead to a better understanding of the modifications necessary to incorporate by Goos-Hänchen and related effects into a ray dynamics around the physically important regime of critical refraction.
In summary, we have shown how the in-out decompositions discussed in this paper lead naturally to expected results for resonant modes of a dielectric circle. We emphasise, however, that the general approach extends to deformed cavities. In particular, we demonstrate that this approach shows promise for a better understanding of effects such as Goos-Hänchen and similarly modified dynamics near critical lines; in the interests of brevity we do not pursue this issue any further in this paper.

Conclusions
We have proposed an approach for separating the boundary solutions of cavity problems into incoming and outgoing components that allows boundary integral equations to be recast as transfer-operator equations. The resulting "shift operator" maps the outgoing component leaving the boundary onto the incoming component arriving back at the boundary having crossed the cavity's interior. The shift operator is completely independent of the boundary conditions and decouples from them the problem of wave propagation through the interior. The boundary conditions are used separately to complete the problem by defining a reflection operator mapping the component arriving at the boundary to the reflected, outgoing component leaving it again towards the interior.
The underlying intuition behind this approach is familiar in the context of semiclassical approximations [1]. We believe, however, that using (8) and (9) rather than ray-based ideas in order to define the in-out decomposition of the boundary solution provides an ideal platform on which to extend these ideas beyond leadingorder approximations. We hope that the exactly defined shift operator will enable the efficient treatment of diffraction and tunnelling effects, for example. The approach also holds promise for more complex problems where small subsystems demanding fully wave-based analysis may coexist with larger structures necessitating ray-based approximations [16,15]. In that context, the ability to perform fully wave-based calculations within a formalism that is suitable also for large-scale ray approximations will allow their efficient integration within hybrid methods.
Equations (8) and (9) lie at the heart of everything calculated in this paper. Their key advantage is in allowing a very natural in-out decomposition to be made while owing nothing to semiclassical approximation. There remains a significant degree of freedom, however, in choosing the singular parts of the Green operators used to construct the outgoing component. The simplest "primitive" choice set out in Sec. 4 appears very natural from the standpoint of ray-based calculations and gives the right leading-order behaviour straightforwardly (see Appendix B). It is singular, however, around criticially incident waves arriving tangentially at the boundary and fails to define appropriately decaying outgoing components in evanescent regimes (see Fig. 7 and the surrounding discussion in Sec. 6.3). We have therefore suggested in Sec. 5 a second approach appropriate to domains with analytic boundaries. Here, we have used asymptotic analysis of the boundary integrals to motivate a decomposition based on integration contours on the complexified boundary. This approach is at once able to provide semiclassical approximation beyond all orders and amenable to exact analysis. It has been shown to work well in providing a solution of model problems with circular geometry in Secs. 6 and 7.
Three outstanding problems remain which deserve further investigation. First, the applications so far carried out have been to the special case of a circular domain for which the shift operator can be fully determined analytically. The regular decomposition is also applicable to more general (analytic) domains, but will require some numerical intervention in such cases to implement fully. Efficient means of performing steps such as the operator inversions in (35) therefore need attention, as do practical approaches to truncation of the complex contours in the presence of chord-length singularities expected for generic boundaries. Second, while we hope that equations (8) and (9) will eventually provide a basis also for the treatment of nonconvex domains, or boundaries with corners, these will respectively require more careful analysis of the nontrivial operator inversion on the right hand side of (36) and the isolation of the singular parts of boundary integrals at corners. Third, there remains some degree of "ugliness" in the treatment of the transition from open to closed modes in the discussion provided for the regularised decomposition. Although the regularised components do not diverge around criticality as do their primitive counterparts, they do still undergo a discontinuous change. This is due to a coalescence and switching of identities of incoming and outgoing components around critical reflection and may therefore be unavoidable, but one would ideally manage the transition in a smoother way.
Acknowledgments SCC acknowledges support from EPSRC under grant number EP/F036574/1 and the hospitality of the MPIPKS, Dresden, where part of this work was completed while participating in the Advanced Study Group Dynamical Tunneling. HBH acknowledges support by the Libyan Ministry of Education. Further support by the EU (FP7 IAPP grant MIDEA) is gratefully acknowledged.

Appendix A. Green operator on a line
The extraction of the singular partĜ sing 0 of the Green operator plays a central role in the decomposition of the boundary functions on which we base the transfer-operator representation of the boundary integral equations. In the case of the primitive decomposition described in Sec. 4, we have definedĜ sing 0 so that it corresponds to replacing the boundary locally by a straight line. This allows the simple, compact representation given in (20), which is derived in this appendix.
LetĜ sing 0 be defined by (19). Before taking the limit ε → 0 in this definition, the operator takes the form of a convolution of the test function µ(s) with a function whose Fourier transform is where the latter form is tabulated in section 6.677 of [45]. Taking the limit ε → 0, the action ofĜ sing 0 on the Fourier transform χ(p) of µ(p) is such that This justifies the formal definition (20) of the singular part of the Green function, along with the integral version given in (23). The branches of the square root obtained for p 2 > 1 also explain the route of the contour in (23) around the branch points p 2 = 1.
We now derive the equivalent expression (24) forĜ sing 1 , beginning with its formal definitionĜ analogous to (19). Note that applying this operator amounts taking a derivative with respect to ε in (19) before taking the limit. Therefore, its action on the Fourier transform ϕ(p) of ψ(s) is obtained by multiplying it by the derivativẽ of (A.1) and then letting ε → 0, after whicĥ which demonstrates (24).

Appendix B. Semiclassical approximation of the primitive decomposition
In this appendix we describe more explicitly some of the semiclassical results claimed in the main text for the primitive approximation. To describe the primitive decomposition semiclassically, we begin with the following rewriting of (25) We next note that applying the operatorP to a propagating boundary function of WKB form with local tangential momentum p = sin α amounts at leading order to multiplication by | cos α|. Here α denotes the angle of incidence of the associated ray, defined with respect to a normal pointing into the domain Ω and we note that real chords on convex billiards satisfy cos α > 0. Away from the diagonal s = s , where the singular contribution fromÎ can be neglected, the Green function can be approximated by For convex cavities, we find then thatŜ 0 can be represented at leading order by the kernel cos α e ikL(s,s ) , (B.1) where α and α are the chord's angles of incidence defined in Fig. B1. One easily finds on using ∂L/∂n = cos α that a corresponding calculation of the kernel S 1 (s, s ) for S 1 gives the same leading form. This justifies (29) for real rays in convex domains. The situation is different, however, for nonconvex domains or for closed or evanescent modes represented by complex rays. Consider first the case of real rays in nonconvex domains. In this case one encounters "ghost orbits" which pass through the exterior of Ω on the way from s to s. In particular, chords which immediately leave s towards the exterior will be such that ∂L/∂n = cos α < 1, in which case the kernels ofŜ 0 andŜ 1 have opposite sign at leading order: Note that ghost orbits can be shown to cancel systematically to all orders in the full secular equation [7], but their presence significantly complicates the present analysis by preventing a simple evaluation of the operator inversion in (28). More generally they prevent us from imposing (35) in a simple way. This is not surprising as we expect the corresponding external problem, discussed at the end of Sec. 5.2, to have nontrivial incoming components in nonconvex problems. We state without proof that complex chords corresponding to evanescent wave propagation beyond criticality also satisfy (B.2). This is reflected also in the vanishing of the right hand side of (50) at leading order when m 2 > (kR) 2 . The result is that ψ − =Ŝψ + is in general smaller than ψ + in the closed region of momentum representation, although, as evidenced by Fig. 7, not exponentially so as the cancellation of the right hand side of (28) only occurs at leading order and nontrivial higher-order contributions remain.
We remark finally that a more explicit comparison with standard representations of the transfer operator can be made if the incoming and outgoing boundary waves are represented by Ψ ± defined in (51) instead of ψ ± and, replacing (25) and (26), correspondingly. The transformed shift operatorsŜ 0 andŜ 1 are then defined bŷ The direct analogue of (28) now yields a shift operatorŜ =P 1/2ŜP −1/2 such that Ψ − =ŜΨ + , and which can be represented by a kernel whose leading semiclassical approximation takes the standard form given in [1].

Appendix C. Alternative treatment of evanescent components
The treatment of evanescent components within the regularised decomposition has been described primarily in this paper while assuming that the endpoint contour component Γ 0 has been chosen to pass to the right of the saddle in Fig. 5(b). This has the advantage of simplifying the presentation of formal results by allowing, for example, condition (35) to hold at leading order semiclassically over both the open and closed subspaces.  Figure C1. The decomposition of contours illustrated here provides an alternative to the convention in Fig 5(b). This is not the only possible choice of Γ 0 , however. As described in Sec. 5, the saddle-point contribution to (31) undergoes a Stokes transition around real values of p 2 > 1 so that it is a priori equally valid for Γ 0 to pass to the other side of the saddle, as shown schematically in Fig C1. Choosing this alternative route provides us with different forms forĜ sing 0,1 andĜ reg 0,1 on the closed subspace, which we describe in this section for the case of the circle: note that, because Γ 1 is now nonempty, the corresponding components ofĜ reg 0,1 are nonzero. We emphasise that, throughout this appendix, the restrictions ofĜ sing 0,1 andĜ reg 0,1 to the open subspace are unaffected by this choice.
We find for this alternative contour decomposition that (this identity is not immediately obvious but follows from manipulation of the forms given above forĜ sing 0,1 andĜ reg 0,1 and use of the Wronskian identities for Bessel functions). In particular, we find from (36) that the shift operator of the interior problem still vanishes identically for closed modes |m| > z: The solution of the interior problem is then not qualitatively changed by adopting this alternative convention and we simply find that, beyond leading order, the construction of the shift operator is complicated somewhat in the closed subspace.
On the other hand, the shift operator for the exterior problem is changed qualitatively. According to (39), the closed components of this operator then satisfy so that, although the external shift operator is exponentially small (J m (z) Y m (x) in the evanescent regime), it does not vanish. A convention in which the shift operator for the exterior problem is exponentially small but nonvanishing may be an attractive way of treating tunnelling effects directly. For convex domains, waves leaving the boundary evanescently towards the exterior may encounter, and be reflected at, caustics some distance from the boundary before returning to it. In the circular case such caustics are encountered at a radius kr = |m| > kR where the outgoing Hankel function undergoes a Stokes transition. In the language of ray asymptotics, although a real ray leaving the boundary never returns, a contributing complex ray may (see [34] for an explicit discussion of such exterior complex rays in the context of tunnelling emission from dielectric cavities). In the convention used in the main text, such returning evanescent waves are simply included in ϕ + . The alternative here allows such returning evanescent waves to contribute to a nonvanishing component ϕ − and to be treated separately.
One final comment is that yet another convention is obtained by definingĜ sing 0,1 as a weighted mean of the singular parts defined in this appendix and those used in the main text. Although space does not allow us to present details here, it is found that (C.1) still holds in such conventions, as do the identities (Ŝ int ) mm = 0 and (Ŝ ext ) mm = (Ŝ 0 ) mm . Furthermore, a symmetrically-weighted mean leads, for example, to operatorsD 0 andD 1 having purely imaginary components in the closed subspace.

Appendix D. Mapping boundary solutions to the interior
In this appendix we summarise how the solution inside the domain can be determined in terms of the in-out boundary solutions. Let extend the Green operators defined by (3) and (4) so that they map boundary functions to functions defined on the interior Ω. Then the interior solution ψ(x) can be written formally ψ(x) =Ĝ 0 µ −Ĝ 1 ψ. Written in terms of the in-out components, this becomes (after some manipulation) If the extended Green operators satisfy the following generalisation of condition (35)   where T (x, s ) −2ikG 0 (x, s ; k) cos α and α is the angle at which the line from s to x leaves the boundary. If we place x on a section parametrised by arclength s and such that the line from s arrives with incidence angle α, this can be approximated more explicitly by the form cos α e ikL(x,s ) , (D.4) generalising (B.1). We note finally that an analogous boundary integral provides the solution ϕ(x) of the exterior problem in terms of the component ϕ + (s) leaving the boundary towards infinity, except that the angles of incidence are defined relative to the outward pointing normal in that case: in the dielectric problem, α is replaced by the angle of refraction β .