Frobenius–Perron eigenstates in deformed microdisk cavities: non-Hermitian physics and asymmetric backscattering in ray dynamics

In optical microdisk cavities with boundary deformations the backscattering between clockwise and counter-clockwise propagating waves is in general asymmetric. The striking consequence of this asymmetry is that these apparently weakly open systems show pronounced non-Hermitian phenomena. The optical modes appear in non-orthogonal pairs, where both modes copropagate in a preferred sense of rotation, i.e. the modes exhibit a finite chirality. Full asymmetry in the backscattering results in a non-Hermitian degeneracy (exceptional point) where the deviation from closed system evolution is strongest. We study the effects of asymmetric backscattering in ray dynamics. For this purpose, we construct a finite approximation of the Frobenius–Perron operator for deformed microdisk cavities, which describes the dynamics of intensities in phase space. Eigenstates of the Frobenius–Perron operator show nice analogies to optical modes: they come in non-orthogonal copropagating pairs and have a finite chirality. We introduce a new cavity system with a smooth asymmetric boundary deformation where we demonstrate our results and we illustrate the main aspects with the help of a simple analytically solvable 1D model.

In microdisk cavities an optical mode is characterized by an azimuthal and radial mode number [24]. The azimuthal mode number m indicates the traveling in clockwise (CW) or counter-clockwise (CCW) direction for m 0 < or m 0 > respectively. Due to the rotational symmetry of the microdisk, traveling waves with m  are degenerate in their complex frequenciesΩ (the imaginary part determines the decay rate). Therefore, one can superimpose them to standing waves in azimuthal direction with either positive or negative parity. Deforming the boundary of a microdisk induces a coupling between CW and CCW components. In a twomode model this can be described with the effective Hamiltonian in the (CCW, CW)-basis [25][26][27][28]  By deforming the microdisk in a way that a mirror-reflection symmetry is preserved, the backscattering coefficientsA and B are equal. This symmetric backscattering implies that the CCW and the CW components of the given eigenvector(3) are of the same size. Both modes are therefore standing waves as illustrated in figure 1(a) for an elliptical cavity.
A deformation that breaks all mirror-reflection symmetries leads in general to asymmetric backscattering A B; ¹ | | | | see [25][26][27][28]. Note that the openness of the system is crucial, otherwise A B * = . Inspecting the structure of the eigenvectors in (3) reveals three interesting consequences: chirality, copropagation, and non-orthogonality. Chirality means that the CCW and the CW components of the given eigenvector (3) are not of the same size. Hence, the optical modes are not standing waves but partially traveling waves [25,26,29,30], a fact that has been confirmed in a recent experiment [31]. Both eigenvectors(3) have the same dominant component and, therefore, both modes have the same preferred sense of rotation, i.e. they copropagate. Finally, the optical mode pair is non-orthogonal since [25][26][27][28] as illustrated for an optical mode pair in figure 1(b). Proposed applications of asymmetric backscattering are an enhanced sensitivity of particle sensors [32] and optical gyroscopes [33], and optical waveguides with special transport properties [34].
For full asymmetry in the backscattering (A = 0 with B 0 ¹ or the other way around), one gets perfect chirality, collinear eigenvectors, and degenerate eigenvalues. These properties indicate an exceptional point (EP). An EP is a point in parameter space of a non-Hermitian system at which two or more eigenvalues and eigenstates coalesce [35][36][37][38]. The physical existence of EPs has been demonstrated by experiments on a number of physical systems as reviewed in [39]. The first direct experimental verification in a wave system has been done for a microwave cavity [40,41]. Later, EPs have been observed in optical microcavities [19,20]. Moreover, EPs play an important role in  -symmetric systems [42]. These systems can be realized by electromagnetic systems which have a balanced arrangement of absorbing and amplifying regions [43][44][45].
In [46] chirality has been predicted in the long-lived intensity distributions of rays inside a spiral-shaped cavity. Moreover, asymmetric backscattering has been directly observed in numerical simulations of the same kind of cavity [25] and later also in other asymmetric cavities [26]. However, it was also shown in [26] that there are cavity geometries where the ray dynamics is not able to feel the asymmetric backscattering; see also [27]. The panels show (a) the intensity of a pair of nearly degenerate orthogonal modes in the symmetric ellipse cavity. The modes are standing waves with stationary mode pattern and without net propagation of light intensities. In (b) the intensity of a nearlydegenerate non-orthogonal mode pair in the asymmetric Fourier-truncated spiral (see (8)) is shown. The modes are traveling waves with stationary mode pattern and net propagation in CW direction, see red arrows. Black arrows illustrate the (C)CW components of the individual modes.
In this paper we show that when the ray dynamics is sensitive to the asymmetric backscattering there is not only chirality in the long-lived intensity distributions of rays but also copropagation and non-orthogonality of certain ray dynamical states. We define these states as eigenstates of the Frobenius-Perron operator  [47] for the time evolution in phase space. We observe that the eigenvalues of  come in pairs and show that at least two eigenstates with large eigenvalues are important for the time evolution. As indications of asymmetric backscattering these two eigenstates are highly non-orthogonal and have a finite chirality meaning they tend to copropagate in a specific (CW or CCW) direction.
This paper is organized as follows. In section 2 we describe the system studied in this paper and fix the notation for phase-space dynamics. Section 3 is divided into five parts where we describe the construction of a finite approximation of the Frobenius-Perron operator (section 3.1), investigate properties of the eigenvalues (section 3.2) and eigenstates (section 3.3) of  , which are then used to construct a two-mode model for asymmetric backscattering (section 3.4) and we compare the results of ray dynamics and wave calculations (section 3.5). In section 4 we present a toy model in 1D which resembles the properties of the eigenstates of the 2D phase-space dynamics. The results are summarized in section 5.

System and phase-space dynamics
Rays inside a microcavity with homogeneous refractive index n propagate on straight lines until boundary reflections that change the direction of the rays. A convenient way to analyze this dynamics is a Poincaré section of phase space with Birkhoff coordinates s p , ( ). Here, s determines the position along the boundary, where a ray with angle of incidence χ is reflected, see figure 2. The coordinate p sin c = is the tangential momentum of the reflected ray at boundary point s. Thus, s p , ( ) define a ray emerging from the boundary. Its time evolution up to the next boundary reflection is determined by a Poincaré map (see figure 3(a) for an example) Since s p , ( ) are canonical conjugated variables, the mapping  is area preserving. In the following, we identify rays with p > < ( )0 as CCW (CW) propagating rays. In addition to the billiard dynamics of the closed system which is described by , the intensity I of a ray is a dynamical variable that decreases at boundary reflections. This loss of intensity is determined by the reflectivity R p ( ), that is a function of the tangential momentum p. For p n 1 < | | and transverse magnetic (TM) polarized light it is given by the reflectivity is unity due to total internal reflection (see figure 3(b)). Thus, the intensity I of the phase-space point s p , ( ) is mapped via to the intensity I˜of the iterated phase space point s p , (˜˜). Since (4) is a discrete version of continuous dynamics, time also has to be considered as a discrete variable that is increasing monotonously with an increment t D depending on the length l D of a ray. Here, the increment t D is a function of initial conditions s p , ( ) in phase space (see figure 3(c) for an example). Thus the time variable is mapped via In the following we focus on a specific smooth boundary given by the radius as function of polar angle f. R determines the size of the system. For numerical studies we use R 1 = . Other parameters can be chosen within a broad range. Here, we fix them to 0.07  = and N 4 = (see boundary in figure 2) so that chaotic ray dynamics is dominant which is illustrated in figure 3(a). The refractive index is set to n 3.0 = leading to partial leakage for p 1 3,1 3 Î -[ ] , see figure 3(b). Since the sum in (8) contains sine terms the system does not contain any mirror-reflection symmetry and r r )for any 0 f . Note that for N  ¥ the non-smooth spiral [15,25,48,49] with a notch of size R  p is recovered. We therefore call the cavity defined by (8) the Fourier-truncated spiral.

Frobenius-Perron operator for microcavities
We discussed the iteration of a point in phase space that represents a single ray in the previous section. Now we focus on the iteration scheme of an ensemble of rays which is given by a density ρ in phase space. The time evolution of ρ is then determined by the Frobenius-Perron operator  as For closed billiard systems  is given by [47] s p s p J s p , , , where the Jacobian J s p , 1 = ( ) for Hamiltonian dynamics.
Since microcavities are open systems that lose intensity through a leaky region in phase space, the Frobenius-Perron operator needs to be adapted. If additionally the non-unique time increment t D is taken into account, the Frobenius-Perron operator becomes [50] s p R p s p , e , , 1 1 where the escape rate κ quantifies the average loss of the intensity in the system. Consequently the Frobenius-Perron operator is sub-unitary with eigenvalues inside the unit circle in the complex plane. The structure of an eigenstate of  is invariant under time evolution but the measure of the eigenstate is reduced after each time step. The long-time dynamics of the system is here dominated by the eigenstate with largest eigenvalue modulus because other eigenstates lose their intensity faster than this one. For billiard systems and symplectic maps the largest eigenvalue turns out to be real and the corresponding eigenstate is positive as observed in [50][51][52][53][54][55][56]. In case of asymmetric microcavities we show that there exist a pair of nearly degenerate eigenstates with large eigenvalues of similar size. We observe that the dynamics is dominated by at least two eigenstates of  .

Construction of a finite approximation of 
Since the Frobenius-Perron operator has infinite dimension one needs proper approximation schemes to get finite and therefore numerically diagonalizable approximations of  . Typically,  is projected onto a suitable set of basis functions which are for instance spherical harmonics as in [52] for an angular momentum map. In a general case it is convenient to divide the phase space into cells and choose characteristic functions of these cells as basis for  . Thus,  is expressed in densities scattered between cells which are calculated by applying the phase-space map. This procedure is known as Ulamʼs method [57] and has been applied (with slight modifications) to the standard map [53][54][55]. In this paper we adapt Ulamʼs method to calculate a finite approximation of the Frobenius-Perron operator for microcavities. Therefore, we calculate not a scattering of density but of intensity between cells with finite reflectivity and include the true-time aspect of the billiard dynamics according to (11).
In the following we explain our procedure in detail: first, the phase space is divided into N grid × N grid cells of equal size labeled with an index. In each cell j random initial conditions s p , (˜). We select the points that are iterated to cell i and calculate the scattered intensity I j i  ( )from cell j to cell i as where N τ is the number of trajectories started in cell j. This procedure is sketched in figure 4.
As shown in figure 3(c) the time to the next boundary collision is quite different for individual initial phasespace conditions. Orbits with large tangential momentum slide along the boundary and need a large number of iterations to get a desired length whereas the length of orbits with small tangential momentum increase fast with the number of iterations. To respect this true-time aspect of the billiard dynamics in the Frobenius-Perron operator we need to weight the scattered intensity I j i  ( )according to the length (or fight time) of the orbits. Inspired by [50] we choose the weighting factor consistent with (11). Here, t ij D is the average time over the N t trajectories for one iteration from cell j to cell i and á ñ · denote an average over cells in one iteration. A simple estimate of the escape rate κ is given by being the measure of the leak and t áD ñis the average time between boundary collisions [50]. Thus the matrix element of the finite approximation of the Frobenius-Perron operator is given by » t trajectories to get statistically satisfying results. We end up with a sparse N grid 2 × N grid 2 matrix ij  as finite approximation of the Frobenius-Perron operator.

Eigenvalues of 
In this section we focus on the properties of the eigenvalues of the Frobenius-Perron operator that we calculated by diagonalizing the finite approximation of  . Due to its sparseness diagonalization is done with Arnoldi method [58] for N 100 grid > . l » -. So the cavity would lose more or less all of its intensity until one can neglect the second eigenstate. Note that almost all calculated eigenvalues even with non-vanishing imaginary part come in nearly degenerate pairs as shown by the magnifications in figure 5.
The two largest eigenvalues of  turn out to be quite stable against variation of N grid as shown in figure 6(a), whereas the convergence of eigenvalues with lower magnitude is slower. Also the difference between the largest eigenvalues seems to saturate for sufficiently large values of N grid (see figure 6(b)).

Long-lived eigenstates of 
In this section we investigate the eigenstates of the Frobenius-Perron operator  . For N 2000 grid = the two eigenstates 1 2 r with largest eigenvalue modulus are shown in figures 7(a) and (b). They are localized on the unstable manifold of the cavity which is the set of points in phase space that converges to the chaotic saddle in backward time evolution. The chaotic saddle itself is the set of points that does not reach the leaky region in phase space neither in forward nor in backward time evolution. For the Fourier-truncated spiral the unstable manifold is shown in figure 7(c). Its computation is done via the algorithm presented in the appendix C of [50] with t 70 * = and t 8 D = . The localization of eigenstates on the unstable manifold is consistent with the energy density obtained by [59] for a symmetric cavity. Note that the chaotic saddle is also an eigenstate of  whose eigenvalue is one but since it is a measure-zero set we are not able to resolve this state with any finite discretization N grid . Note that the second eigenstate shown in figure 7(b) is negative in the upper half of the phase space. Therefore it cannot be interpreted as ordinary intensity. But in combination with the first eigenstate it is possible to construct linear combinations from 1 r and 2 r which are positive and remain positive under time evolution. Since the cavity presented in this paper does not contain any mirror-reflection symmetry, the weight   . This indicates that both eigenstates mainly propagate in CW direction. For a pair of optical modes the chirality α is connected to their overlap S by S S 2 1 a = + ( ) [26]. For Frobenius-Perron eigenstates this equation is still valid with S , 3.4. Two-mode approximation from  In the following we construct a 2×2 model for the phase-space dynamics of (C)CW regions. Therefore, proper basis states representing pure (C)CW propagation in phase space are needed. Wave calculations [26] have shown that these basis states are given by linear combinations of optical modes with nearly degenerate complex frequencies.
Since the Frobenius-Perron operator for ray dynamics has two relevant eigenstates, we can construct dominantly propagating states in the similar way by , . ) as e.g. in the limaçon cavity [12] there exists a small splitting of largest eigenvalues but their eigenstates are (anti-)symmetric with vanishing chirality (18). .
The insets show the corresponding states in phase space from low intensity (white) to high intensity (red/gray).
In the case of full asymmetry in the intensity backscattering A where both eigenvalues and eigenvectors are degenerate. But as already seen in the full numerics (see figure 5) we have nearly degenerate eigenvalues. Therefore, we get A B , , 0.925 249, 5.4 10 , 7.9 10 The off-diagonal elements determine the backscattering of intensity from CCW to CW and vice versa. In our example B A  indicate that the backscattering from CCW to CW is much larger than in the reverse process. The eigenstates of the Hamiltonian H are as in (3). Their normalized overlap

Comparison to wave calculations
In this section we compare our results of ray dynamics with wave calculations. Therefore, we use the boundary element method [60] to solve the 2D mode equation for TM boundary conditions. Here, the function kR Y and its normal derivative are continuous at the boundary of the cavity. As shown in figure 1(b) we typically get pairs of modes with nearly-degenerate values of kR. The phase-space representation of such a mode kR Y giving via its Husimi function s p , kR  ( ) [61], is comparable to eigenstates of the Frobenius-Perron operator. For closed billiard systems kR is real and the average over all Husimi functions up to a given kR converges to the invariant measure for kR  ¥ [62]. In contrast microcavities are open systems with complex kR. In [63,64] it was shown that an average around a large RekR over many low-loss modes should correspond to the ray dynamics. Therefore, we randomly determine N 183 kR = optical modes with RekR ä [70, 72] and quality factor Q = −RekR/(2ImkR) > 2000. From these modes we calculate an averaged Husimi function  via As shown in figure 9 the averaged Husimi function gives a nice correspondence to Frobenius-Perron eigenstates 1 2 r . The chirality 0.952 311  a = of the averaged Husimi function is calculated according to (18) and fits quite nice to the values obtained from 1 2 r . Note that optical modes with RekR≈70 are not able to resolve the very thin structures of the unstable manifold or the eigenstates of  , since an associated effective Planck constant h=2π/(nRekR) is too large compared to the resolution N 1 grid 2 µ of the finite approximation of Frobenius-Perron eigenstates.

Asymmetric 1D model
In this section we construct a toy model which reduces the complexity of the dynamics in a microcavity to a simple 1D system. Therefore this 1D-toy model is a minimal example of a dynamical system with asymmetric backscattering between two weakly coupled regions. Furthermore the eigenvalue problem of the Frobenius-Perron operator is solved analytically in this case and it is shown that the toy-model eigenstates have quite similar properties as the eigenstates discussed in section 3.3 including non-orthogonality, chirality, and copropagation. ] with a function R x ( ) (see figure 10(c)). In contrast to the reflectivity from (5) we here assume the function R x ( ) to be asymmetric in the region of exchange a a , Here A and B quantify the intensity which is scattered from 1, 0 -[ ) to 0, 1 [ ]and vice versa. The intensity which remains in either a 1, Thus, the dynamics of a state ρ for one time step is determined by a a for 1, 0, , for , 0 , 1 , 28 We solve the eigenvalue problem r r = W  with the following ansatz x a a for 1, 0, , for , 0 , 1 , 29 where C 1 and C 2 are constants to be determined. Applying  given by (28) and compare both sides of the eigenvalue problem yields This equation is the equivalent to the eigenvalue problem corresponding to the two-mode Hamiltonian in (1). Thus, the Frobenius-Perron operator has two eigenvalues as in (2) and two eigenstates 1,2 r with coefficients C C , 1 2 ( ) given by (3). In the following example we specify x a a x a 0.  Linear combinations of the eigenstates according to (19) lead to states that are dominantly localized either on the left or right side of the exchange interval as shown in figures 12(a) and (b). Therefore, they are identified with CW and CCW propagating states respectively. To confirm the asymmetry of the backscattering we simulate the time evolution of 10 7 points with an initial intensity distribution either in CW or CCW region. Thus, we get intensities which are scattered from 1, 0 -[ ]to 0, 1 [ ]and vice versa. As shown in figure 12(c) more intensity is scattered from CCW to CW region than in the reverse process, while the loss of intensity from the initial region is almost equal for both CW and CCW.

Summary
We have demonstrated the construction of a finite approximation of the Frobenius-Perron operator  for optical microcavities. A new class of boundary deformations, the Fourier-truncated spiral, which destroy mirror-reflection symmetry was introduced. Due to the asymmetry of this system the largest eigenvalues of the Frobenius-Perron operator are nearly degenerate. The two corresponding eigenstates are important for the time evolution of intensities inside the cavity. The eigenstates of  reveal non-Hermitian properties of the ray dynamics which are known from optical modes of asymmetric microcavities: they are pairwise non-orthogonal and they are favored to copropagate in CCW direction, which is measured by a finite chirality.
We compared our results from ray dynamics to results from wave calculations. A nice agreement between eigenstates of the Frobenius-Perron operator and an averaged Husimi function of low-loss modes is found.  In addition, we presented a 1D toy model where the Frobenius-Perron operator is diagonalizable analytically and the resulting eigenstates have similar properties to the eigenstates of the Frobenius-Perron operator of the microcavity.
Note that for the appearance of non-orthogonality, chirality, and copropagation of long-lived Frobenius-Perron eigenstates it is necessary that the microcavity is treated as an open system with a leaky region in phase space and that the asymmetric backscattering is visible for the ray dynamics. This implies that the boundary curve has a certain asymmetry which leads to mixing/chaotic ray dynamics and a coupling of CW and CCW regions in phase space. Counterexamples are not simply connected cavities like a microdisk with external scatterers [27] or billiards of constant width [26,65]. In these cases, asymmetric backscattering occurs wave mechanically (e.g. by dynamical tunneling) but not in the ray dynamics. But for simply connected and asymmetric microcavities the non-orthogonality, chirality, and copropagation of long-lived Frobenius-Perron eigenstates should be a general feature and we observed it also in the asymmetric limaçon, the spiral and asymmetrically connected half circles and also for TE polarization and lower refractive indices. ]. This sub-matrix has ones in the counterdiagonal. Finally, the matrix R determines the partial leakage of the system and has only diagonal elements R x R ii i = ( ). Because of this matrix R the representation of the Frobenius-Perron operator F becomes sub-unitary with eigenvalues inside the unit circle in the complex plane.