Perturbation approximation for higher modes in nearly regular two-dimensional cavities

A perturbation theory for weakly distorted regular cavity which has classical ray trajectories lying on invariant tori, is constructed to a higher perturbation order, than for the general case. This is possible because of a special structure of semi-classical eigenvalues for integrable Hamiltonians. The perturbation magnitude here has an order of a characteristic wavelength of a mode instead of usual wavelength square. The results are expressed in solutions of the Hill equation. The set includes modes localized along stable periodic ray trajectories; scar modes, corresponding to unstable periodic trajectories; weakly distorted modes of regular cavity, and intermediate cases. The application of the method to square, circular and elliptical cavities is outlined.


Introduction
Finding the modal structure of an optical cavity is important both for applications and from the purely theoretical viewpoint. Recent interest in the topic is related to the theory of microresonators (Cao & Wiersig, 2015;Harayama & Shinohara, 2011). The question is also relevant for quantum chaos manifestation in dynamical billiards (Reichl, 2004). Mathematically, the modes of two-dimensional (2D) cavity are eigenfunctions of Laplace operator in a 2D domain with appropriate boundary condition. For higher modes, numerical calculations are slow and give little insight. The exact analytical solutions are known only for a small number of shapes: circular, elliptic and rectangular among them.

Conformal mapping and perturbation operator
The usual approach to obtaining a perturbation matrix for slightly deformed resonator is based on the Hadamard formula for Green function of a perturbed domain (Bruno & Reitich, 2001;Dubertrand et al., 2008). For two dimensions, another possible method is the conformal mapping (Leonhardt, 2006;Reck, Thomsen, & Hansen, 2011;Schinzinger, 2003). If an analytic function of complex variable, w(z), maps a simple shape to a cavity boundary, the solution of Helmholtz equation in the resonator with complicated boundary is mapped to the solution for the simple resonator filled with a medium which has an effective refractive index distribution |w � (z)| 2 .
Consider an analytic function: w(z) = u(z) + iv(z) of z = x + iy, with values in the region having a boundary Γ. Consider a real function: f (u, v), which is zero on Γ (Dirichlet boundary condition) and obeys the Helmholtz equation: with a Laplace operator Δ u,v = 2 u + 2 v . If Γ is an image of a regular domain boundary S upon w(z), then f (x, y) = f (w(z)) is a function on a regular domain, which has zero value on a boundary S, and the Equation (1) in new variables becomes: The approximation we are interested in is obtained for a small distortion of a shape, i.e. the function w(z), which is with a small addition, |d(z)| ≪ z, |d � (z)| ≪ 1. If we consider eigenfunctions having eigenvalues close to a large number 0 , such as | 0 − | ≪ 0 , we can expand the right-hand side of Equation −Δ x,y f (w(x, y)) = f |w � (z)| 2 .
(3) w(z) = z + d(z), (2), to the first order in d ′ and 0 − , which gives a self-adjoint Schrödinger-type spectral problem on a regular domain: In this case, the unperturbed operator is the Laplace operator A 0 f = −Δf with known eigenfunctions f 0 n,p and eigenvalues 0 n,p . The perturbation operator P is: Simple perturbation methods work for perturbation strength smaller than a distance between levels, which,using the Weyl's theorem on level spacing, gives a usual 2 0 d∕L ∼ 4 ∕L 2 , or d ∼ K −2 ∕L estimation of the Introduction.
However, as can be seen below, simple domains, for which analytic solutions exist, have special spectral structures, corresponding for a semi-classical regime to periodic ray trajectories, and in this case the perturbation problem can be approximately solved in special functions for d ∼ K −1 .

Perturbation for regular cavities
The Laplace operator can be regarded as a quantization of Hamiltonian task of a classical movement of free particle in a billiard, which has a shape of the cavity (Reichl, 2004). If the classical Hamiltonian problem is integrable, action-angle variables exist, and the Hamiltonian is a function of two actions H(I 1 , I 2 ). For the semi-classical approximation, energies (mode eigenvalues) in this case are obtained by setting I 1,2 = q 1,2 + 1 4 1,2 , where q 1,2 are integer mode numbers, and 1,2 are fixed integers, depending on the properties of the trajectory (Maslov indices) (Berry & Tabor, 1976;Waalkens, Wiersig, & Dullin, 1997). These points locally form a square lattice in the plane I 1 , I 2 . The typical case is shown in Figure 1.
Let us fix a line in the lattice, parallel to (m, −l) vector. Consider the situation when it is tangent to the Hamiltonian level line H(I 1 , I 2 ) = const, and m, l are small co-prime integers. Fix local coordinates n 1 , n 2 with the origin in a lattice point close to the contact point. Now introduce another set of coordinates in a plane: This transformation is equivalent to a rotation and changing a scale by √ m 2 + l 2 factor.
For big mode numbers q 1 , q 2 , the curvature of Hamiltonian level lines becomes small, and the energy in new coordinates is expressed to the lowest order in s 1,2 as: where the relation 2 ≫ 1 holds. This gives characteristic arcs in dependence of eigenvalues on the mode number n 1 . The condition that tangent to a level line is parallel to a (m, −l) vector means that corresponding classical trajectory is periodic (Berry & Tabor, 1976). Numerical illustration for a circle is given below (see Figure 4).
Since in the semi-classical region, the energy difference between consecutive arcs is much bigger than the typical energy difference within a single arc, the perturbation, if it is small enough, can be limited to one arc. The energy difference between arcs scales as √ , and thus the perturbation magnitude is limited here by d ∼ K −1 .
If f q1,q2 are eigenfunctions of Laplace operator in a regular domain, we fix a set of modes with the same s 2 : The Laplace operator limited to this set is diagonal: and = H + s2 is the phase, which depends on s 2 value. It can be taken 0 ⩽ < 1 by shifting the origin. The perturbation matrix elements between b k 1 , b k 2 for semi-classical regime depend to the first approximation on the difference |k 1 − k 2 | only. Thus, when k 1 , k 2 are much smaller, than q 1 , q 2 , This gives an approximately diagonal-constant symmetric perturbation matrix. The eigenfunction is where k are components of eigenvector for the infinite matrix given by Equations (9) and (10). Consider a generating function expressed by the Fourier series with the eigenvector components: The action of a matrix on k is equivalent to the action on G( ) of the operator where (7) (s 1 , s 2 ) = 0 + 1 (s 1 − H ) 2 + 2 s 2 , is a 2 -periodic function. Together with the boundary condition G( + 2 ) = G( ) exp(−i2 ), this gives a spectral problem for Hill equation.
In the simplest case, which is easily realized for square and circle, only P 1 ≠ 0, and the Hill equation is reduced to Mathieu equation: where d = P 1 ∕ , and is a rescaled eigenvalue: = ( − 0 )∕ .
For big | d |, the equation can be studied in WKB approximation. Let us consider first = 0, − < < range, and negative d. If d > 0 the solution is obtained from the negative d case by taking − instead of , and the spectrum remains unchanged. Then, for a symmetric function: and for antisymmetric one: These equations are valid between classical turning points, i.e when the expression under the square root is positive. When it is negative, the eigenfunction exponentially diminishes into classically prohibited region of .
The lowest eigenvalue is approximately ≃ −|2d|. For low eigenvalues, the cosine function around the minimum can be approximated by a parabola. This gives a quantum oscillator problem, which has a solution in Hermite functions. The Fourier transform of Hermite function is also Hermite function, which gives k values. If ≫ |2d|, the solution is approximated by slightly perturbed sine/ cosine function, and the corresponding eigenvalue is proportional to a square of a mode number.
Note, that for low eigenvalues and = 0, consequent symmetric and antisymmetric eigenvectors have equally spaced eigenvalues, but for ≫ |2d| pairs are nearly degenerate. The level ≃ |2d| separates these two regimes. Around this value, the distance between consequent eigenvalues diminishes. Also note that the WKB approximation solution close to turning points behaves as The approximate density of modes for eigenvalues around can be obtained using WKB method, and it is expressed in elliptic integrals: . For the case of > 2|d|, the same procedure gives with q + = 1∕q − . The distance between consecutive eigenvalues in function of eigenvalue itself is then given by These equations describe quite well the numerical results (Figure 2(b)).
The behavior of eigenvalues and eigenvectors k is shown in Figures 2 and 3. For these plots, the perturbation matrix of Equations (9) and (10) was taken, and its spectrum and eigenvectors were found numerically. The physical meaning of the resonator mode corresponding to the solution is more easily understood for a circular resonator, and it will be considered below.
Smooth potentials P( ), where P q coefficients rapidly diminish in function of q will result in qualitatively similar behavior which can be investigated by WKB approximation. However, some distortions of boundary produce potentials which are not smooth. One case is when the right angle in the rectangle changes, or when the curvature of a smooth boundary becomes infinite under mappingthe conformal map then has a singularity.
Another problem exists if the trajectories have different topological properties, as it is observed for the ellipse (Sieber, 1997;Waalkens et al., 1997). In this case, the phase space is divided into regions, where Maslov indices 1,2 are different. The lattice at the boundary of the regions is not regular in this case, and the matrix of a perturbation that mixes states from both sides of the dividing line in the action space cannot be constructed with the outlined approach. It still works however in the regions of action plane far enough from such dividing lines.

Square cavity
For a square, calculations of the previous two sections can be done explicitly. Take 0 < x ⩽ 1, 0 < y ⩽ 1. The Hamiltonian (if mass is 1/2) is H = p 2 x + p 2 y , with standard momenta. In function of two actions it is expressed as (Itin & Neishtadt, 2003): and the actions are proportional to absolute values of corresponding momenta I x,y = |p x,y |∕ . The level lines are thus circles, the curvature is constant. The gradient vector is parallel to (|p x |, |p y |), and the condition for periodic trajectory is |p x |m = |p y |l.
The conformal maps can be chosen as: where all c m are real, and |d � (z)| ≪ 1 inside the square. Then the square distortion under the conformal map is such that the sides x = 0, x = 1, y = 0 remain the same, but y = 1 side is transformed in a limit of small perturbation into the curve with equation If only c 1 ≠ 0, this is equivalent to one half of the ripple billiard (Li, Reichl, & Wu, 2002). If a small smooth boundary perturbation y(x), which has y � (0) = y � (1) = 0 is given, it can always be expanded in a series Equation (9).
The matrix elements of perturbation from the Equation (4), which correspond to only one term in Equation (25) are then given by the integral over unit square They are equal to if n 1 = n 2 ± m, n 1,2 > m, and q = p 1 − p 2 . If n 1 ≠ n 2 ± m, the matrix element is zero.
Let us choose the mode which corresponds to a periodic trajectory with short period, such as n 0 , p 0 are related as n 0 m = p 0 l (m, l are small co-prime integers). Then select a subset of modes b k , along s 1 axis, which passes through n 0 , p 0 point. It is given by Equation (8) with q 1 = n 0 , q 2 = p 0 .
(s 1 , s 2 ) = 0 + 1 (s 2 1 + s 2 2 ) + 2 s 2 , This corresponds to considering only one arc. It can be done, while c m ≲ 1∕ √ 0 ∼ K −1 , which means that perturbation is smaller, than a characteristic distance to the next arc. The diagonal matrix elements of Equation (4), limited to b k are: If we have only one term in the sum of Equation (25) c m ≠ 0, the only non-zero off-diagonal elements are: If p 0 ≫ l(2k + 1) this is reduced to Thus, the problem is reduced to the solution of Mathieu equation of previous section. The same procedure with simple modifications can be used for a rectangular cavity 0 < x < L x , 0 < y < L y .

Circular cavity
Consider a unit circle with cylindrical coordinates r, , such as x = r cos( ) and y = r sin( ), For a circle, actions are given by (Ree & Reichl, 1999): where m = 1∕2, and p , E are the angular momentum and the energy, which are conserved.
The dependence I r (I ) with H = E = const, as expressed by Equation (33) is a level line. It is seen, that I r ∕ √ E is an energy-independent function of I ∕ √ E , thus all level lines are obtained by scaling any of them in two directions. The curvature does not change sign, and tends to infinity for I → 1. The condition for periodic trajectory is that a derivative dI r ∕dI in Equation (33) is a rational number, which gives: The Laplace operator eigenfunctions for a unit disk with Dirichlet boundary condition are: where J n are Bessel functions, and a n,p their zeroes. Normalizing constants n,p are calculated with scalar products of Bessel functions by themselves. Corresponding Laplace operator eigenvalues are 0 n,p = a 2 n,p .
To obtain conditions for periodic trajectory in terms of Bessel functions, it is necessary to place in Equation (34) eigenvalues of corresponding operators, which are n for angular momentum, and a n,p for square root of energy. Thus the periodic trajectory condition is: (35) f n,p (r, ) = n,p J n (a n,p r) exp(in ), n∕a n,p ≈ cos( l∕m).
This condition also can be derived directly from well known asymptotic of large order Bessel function (Watson, 1966): The Bessel function zeroes close to 1,000 are shown in Figure 4. The picture clearly demonstrates arcs similar to those obtained for a case of rectangle. The equations for arc parameters can be obtained either from Hamiltonian, by calculating curvatures, or directly from the Equation (37).
The relevant conformal mapping is: The perturbation for one term with real c m is a multiplication by −2 0 (m + 1)c m cos(m )r m .
It is easy to see, that to the first order in |c m | the boundary Γ is given by Thus, by performing Fourier transform one can reconstruct c m coefficients if the shape of weakly distorted circle is known.
The perturbation matrix for only one term c m in the vicinity of resonance has the same three diagonal structure, as for the case of a square.   with 12 = n1,p1 n2,p2 . The angular part of this integral differs from zero only if n1 = n2 ± m, as it was for the square, but the radial part is not expressed in elementary functions. However, if |p1 − p2| ≪ p1, p2, these integrals do not depend strongly on p1 + p2, as for the rectangular case.
Thus, For r → 1, all Bessel functions involved are close to cosine functions with similar periods and amplitudes, and they have zeroes for r = 1. The sign of them can be the same, or alternate with k depending on l value in an arc. Thus, the angular intensity distribution for a mode g (r, ) close to the circumference r → 1 nearly repeats the modulus square of G function from Equation (18)

Elliptical cavity
The case of ellipse is more complicated, and we consider it only briefly. Elliptic coordinates naturally appear within the conformal mapping framework. Consider a function defined on the rectangle 0 < x ⩽ L x , 0 < y ⩽ 2 . Its real and imaginary parts are: which gives usual expressions for elliptic coordinates. The Equation (42) maps the rectangle onto the inner part of the ellipse, the side x = 0 is projected onto the cut between foci. The Laplace operator Equation (2) in these coordinates becomes and it is separable, f (x, y) = F x (x)F y (y). The equations for angular part F x (x) and radial part F y (y) are Mathieu and modified Mathieu equations, where = 2q, and is separation constant. The solution of Mathieu equation has to be 2 periodic; for a given q this gives a set of possible values. The eigenvalue is fixed with a boundary condition for F x and continuity across a cut (see Sieber, 1997, Waalkens et al., 1997 for details).
A simple conformal map can be constructed in a form of Taylor series of Equation (38). General theorems ensure that a smooth boundary perturbation produces a rapidly converging series. However, calculation of perturbation matrix elements in this case will require 2D integrals over the ellipse.
−Δf = (cosh(2x) − cos(2y))f , It is possible to reduce the problem of finding an appropriate conformal mapping to the maps, which slightly deform a rectangle, (z). The ellipse distortion can be then obtained by applying the transformation of Equation (42). Since x variable is cyclic, terms exp(mz) are involved. The map, that does not produce a discontinuity in |w ′ | 2 across a cut has a form. and c m are complex coefficients. Consider, that only c m ≠ 0. For real c m , the x = L x side is displaced in x-direction according to c m sinh(mL x ) cos(my). For purely imaginary c m = id, the displacement is −d cosh(mL x ) sin(my). Note that in this case, all four sides of the rectangle are distorted, and the cut between the foci is no more straight. However, the derivative absolute value is continuous across a cut.
It is seen, that if the boundary displacement in elliptic coordinates is known, x(y) = L x + d(y), then coefficients c m can be determined by Fourier transform, similar to the case of a circle.
The calculation of perturbation matrix elements is then reduced to the estimation of four 1D integrals.
The Hamiltonian for elliptical billiard, action variables, and correspondence of classical and quantum mechanical solutions are discussed in detail in Waalkens et al. (1997), Sieber (1997. The important difference with a circular case is that for elliptical billiard two types of movement are possible, and in the action plane I 1 , I 2 two regions are seen, divided by a straight line I 2 = I 1 . The transition across the line is not continuous. For a correct semi-classical approximation the square lattice in the vicinity of a transition line become distorted as well (Sieber, 1997). Thus, the discussed approximation cannot be built in the vicinity of the classic trajectories passing close to the foci, and special consideration is required there. These trajectories under perturbation form a homoclinic tangle in a phase space. Numerics demonstrate, that if the states far from dividing line are considered, the matrix elements diminish rapidly with |k 1 − k 2 |, and the spectrum becomes wider if we approach the dividing line. More detailed investigation is out of the scope of this paper.

Relation to numerical solution
To illustrate the application of the method, we show here the eigenvalue calculation for the slightly distorted unity circle. We study the arc, which has a maximal eigenvalue equal to 10,624.00 for n = 50, corresponding to m = 3. The conformal mapping function is w(z) = z + 0.01z 4 . The numerical calculation of eigenmodes and eigenvalues for distorted circle is made with a freely distributed "FreeFem++" package based on the finite element method, with a mesh of approximately 200,000 triangles. To build the perturbation matrix, 14 circle eigenmodes are taken around the maximal one, their eignvalues are given by the Bessel function roots. The off-diagonal elements of the Equation (40) are estimated at the central eigenvalue. For K ≈ 100, the semi-classical approximation is not very good, in particular the perturbation matrix coefficient vary 20 percent from one edge of the arc to another. We could not work with higher K values because of finite element software limitations. However, the sequence of modes similar to those of the Figure 5 is easily identified, except for a scar mode Figure 5(d), which is not seen because of a small K number. We compare eigenvalues for 10 higher modes. The first Hermite-Gauss mode (Figure 6(a)) has a calculated eigenvalue 10, 807 ± 2, the approximation gives 10,804. For the fourth mode (Figure 6(b)) we obtain 10,603 and 10,609, the sixth mode (Figure 6(c)) has 10,448 and 10,451, and the tenth (Figure 6(d)) has 10,322 and 10,324 values, respectively. For the rest of modes, the results are similar. The error of the finite element method was estimated by calculating eigenmodes of the circle. The error is due to a finite mesh, but making a finer mesh resulted in program instability. It is seen, that the method gives a reasonable approximation even for the case not very close to a semi-classical one.

Discussion and conclusions
The results have close relationship to the initial stages of KAM chaos development in a slightly distorted regular billiard. The wave structures corresponding to initial stages of KAM chaos in ray optics are obtained. In fact, the reported construction seems natural for investigation of classical behavior of trajectories, and their characteristic features, such as splitting. We do not consider these questions here in any detail, because this has no direct relation to the topic. However, keeping in mind this correspondence is useful for predicting the character of solutions. The relation to KAM chaos makes estimations of validity range for approximations a laborious task. Note that typical perturbations include, to some extent, higher harmonics, thus, modes corresponding to trajectories with big periods can become distorted in complicated ways. Higher perturbation terms also become quite involved. The validity range probably has to be investigated independently for each particular case.
The advantage of the approach is that it gives, in a unified way, different types of modes observed in cavities: nearly regular ones, semi-classical modes around stable trajectories, and scar modes. The last type is especially difficult for analytic methods. The approach also seems to be potentially useful for modes corresponding to a classical KAM homoclinic tangle, which is a question not completely understood till now.
The main purpose of this paper is introducing a perturbation method, which works to a higher order, than usual ones. The examples of the paper, from our point of view, demonstrate its usefulness in practical calculations, as well as for understanding the general character of solutions.